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ABSTRACT 

Two  dereverberation  techniques  are  applied  to  synthetic 
seismic  data  and  their  performance  in  removing  water  column 
multiples  is  evaluated  and  compared.   One  method  is  an  applica- 
tion of  homomorphic  deconvolution  and  the  other  utilizes 
linear  estimation  based  on  a  minimum  mean  square  error  criterion 
The  analytical  formulations  of  both  methods  are  discussed. 
Performance  is  evaluated  in  terms  of  three  criteria:  percent 
of  multiple  energy  removed,  percent  of  signal  (reflector) 
distortion,  and  visual  improvement  of  the  data.   Results  are 
presented  which  represent  the  performance  of  both  algorithms 
for  a  range  of  environmental  and  signal  processing  parameters 
including  white  noise  level,  multiple  coherence,  reflector/ 
multiple  overlap,  filter  parameters  and  water  column  travel 
time  estimate.   The  techniques  are  found  to  have  comparable 
effectiveness  on  the  synthetic  data;  however,  indications  are 
that  homomorphic  dereverberation  has  greater  potential  in 
shallow  water  applications  while  the  linear  technique  appears 
to  be  more  efficient  for  deep  water  data. 
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CKAPTER  I 
INTRODUCTION  AND  PROBLEM  STATEMENT 

The  geologic  structure  of  the  earth  beneath  the  seafloor 
is  most  often  determined  by  seismic  profiling.   The  procedure 
generally  involves  excitation  of  an  impulsive  acoustic  source 
near  the  sea  surface  and  recording  of  the  reflected  earth 
response  with  a  hydrophone  array.   Low  frequency  sound 
penetrates  the  bottom  and  propagates  in  the  substrata  with 
reflections  occurring  at  discontinuities  in  the  acoustic 
impedance  of  the  earth.   The  thickness  and  density  of  sub- 
bottom  layers  may  be  estimated  from  the  reflected  acoustic 
signal,  or  seismogram. 

The  earth  is  modelled  as  a  discrete  layered  medium  with 
distinct  interfaces  for  most  seismic  applications.   This 
assumed  structure,  while  not  strictly  accurate,  has  led  to 
good  processing  results  in  practical  seismic  work  and  has  the 
additional  advantage  of  being  analytically  tractable.   We  shall 
employ  this  assumption  throughout  the  present  analysis.   A 
detailed  description  of  the  earth  model  used  is  given  in 
Chapter  III. 

The  amount  of  energy  reflected  at  a  discontinuity  is 
ideally  measured  by  the  reflection  coefficient, 


r2C2    rlCl 
r2c2  +  rlCl 
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where  r,  ,  r„  and  c,,  c2  are  the  densities  and  acoustic 

velocities  in  media  1  and  2 .   Here  the  sound  propagates  from 
medium  1  to  medium  2  and  normal  incidence  has  been  assumed. 
The  assumption  of  normally  incident  plane  waves  is  generally 
made  in  single  channel  (one  receiver)  seismology  since  it 
leads  to  analytical  simplicity  and  appears  to  be  reasonably 
accurate.   In  utilizing  this  simplified  model  we. have  ignored 
near  field  effects  and  spreading  losses.   These  are  not  of 
major  importance  in  single  channel  dereverberation  and  their 
inclusion  would  unnecessarily  complicate  the  earth  model. 

It  is  evident  from  the  reflection  coefficient  expression 
that  large  reflections  occur  at  points  of  significant  change 
in  impedance.   The  largest  impedance  discontinuity  encountered 
by  seismic  signals  occurs  at  the  water-air  interface  where  R 
is  nearly  -1.   The  water-bottom  interface  is  also  a  strong 
reflector  in  most  cases.   Thus,  the  water  column  becomes  a 
reverberating  channel  wherein  a  significant  portion  of  the 
source  energy  is  trapped.   Repeated  reflections  from  the  bottom, 
or  multiples,  are  received  at  intervals  corresponding  to  the 
two-way  travel  time  of  sound  in  the  water  column.   Deeper 
reflections  from  the  substrata  are  masked  by  multiples  when 
their  arrivals  are  nearly  coincident  in  time.   Since  propaga- 
tion losses  are  much  smaller  in  water,  the  energy  in  the 
multiples  is  usually  large   compared  to  that  in  the  deeper 
reflectors.   Thus,  water  column  multiples  form  an  unwanted 


-7- 


component  of  the  seismogram. 

The  reflection  coefficient  expression  indicates  that 
all  earth  layers  also  introduce  reverberation.   Except  in  some 
shallow  water  situations,  internal  multiples  are  generally 
not  a  serious  problem  for  two  major  reasons.   Most  importantly 
the  acoustic  attenuation  in  the  earth  is  much  greater  than 
that  in  water  so  that  very  little  internal  multiple  energy 
is  actually  returned  to  the  receiver.   Secondly,  the  reflection 
coefficients  at  earth  layer  boundaries  are  usually  small 
compared  to  those  at  the  surface  and  seafloor  so  that  a 
relatively  small  part  of  the  incident  energy  is  actually 
trapped. 

Figure  1  shows  a  synthetic  seismogram  with  strong 
multiples.   Response  amplitude  is  measured  on  the  ordinate 
and  travel  time  in  seconds  on  the  abscissa.   The  first  large 
signal  component  is  the  bottom  reflection  at  one  second  of 
travel  time  or  about  750  meters  water  depth.   The  first 
multiple  is  an  attenuated,  phase-inverted  replica  of  this 
reflection  at  two  seconds.   Note  that  the  return  from  a 
reflection  horizon  at  two  seconds  travel  time  would  coincide 
with  this  multiple  and  be  obscured.   The  overall  periodicity 
of  the  multiples  is  apparent  in  this  plot.   Actual  reflectors 
occur  at  1.5,  2.2  and  2.9  seconds.   Figure  2  shows  the  seismic 
environment  which  would  produce  such  a  seismogram. 

The  first  practical  multiple  analysis  and  dereverberation 
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Figure  2   Earth  structure  v;hich  would  produce  seisraogram 
like  that  shown  in  Figure  1. 
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algorithm  was  proposed  by  Backus  [1],   His  approach  was  to 
characterize  the  water  column  as  a  sharp,  ringing  filter  with 
an  impulse  response  composed  of  a  weighted  sum  of  delayed 
impulses*   The  weights  and  delays  are  determined  by  the  bottom 
reflection  coefficient  and  the  water  column  travel  time 
respectively.   This  model  leads  to  a  three-point  operator  with 
elements  spaced  at  intervals  corresponding  to  the  two-way 
water  travel  time.   Implementation  of  the  Backus  filter  requires 
estimation  of  the  bottom  reflection  coefficient  and  water 
column  travel  time.   Several  aspects  of  the  performance  of 
this  method  are  discussed  in  Chapter  II. 

Spatial  processing  has  also  been  used  to  reduce 
multiples  [2].   Spatial  schemes  normally  require  multichannel 
arrays  of  large  physical  extent  which  can  be  effectively  focused 
to  discriminate  against  reverberation.   Such  systems  are  widely 
used  and  quite  effective  in  shallow  water  but  their  costs,  both 
for  hardware  and  data  processing,  are  very  high.   Hence,  there 
is  still  a  need  for  time  domain  multiple  removal  techniques  in 
deep  water  situations  and  for  single  channel  systems. 

Two  techniques  have  recently  been  applied  to  seismic 
multiple  removal  with  demonstrated  success.   The  first,  an 
inverse  filter  algorithm  based  on  a  tapped  delay  line  model,  is 
due  to  Baggeroer  [3],  and  is  referred  to  hereafter  as  the  TDL 
filter.   A  tapped  delay  line  is  simply  a  realization  of  the 
time  domain  convolution  of  a  signal  and  a  gapped  operator  [4]. 
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Results  of  this  method  have  proven  superior  to  those  of  three- 
point  filtering,  at  least  in  some  deep  water  situations. 

A  nonlinear  filtering  scheme  using  homomorphic  deconvolu- 
tion  has  been  applied  to  several  aspects  of  seismic  signal 
processing.   The  use  of  this  method  for  dereverberation  has 
been  demonstrated  by  Stoffa,  Buhl  and  Bryan  [5].  The  homomorphic 
transformation  is  essentially  a  mapping  from  convolution  to 
addition  so  that,  after  transforming,  deconvolution  can  be 
accomplished  by  simple  linear  filtering.   Seismic  dereverbera- 
tion appears  theoretically  to  be  a  very  promising  application 
of  homomorphic  deconvolution  because  of  the  distinct  properties 
of  seismic  signal  components.   The  method  has  not,  however, 
been  fully  evaluated  or  widely  used  in  practice. 

The  motivation  for  this  study  arises  from  the  disparate 
theoretical  mechanisms  by  which  these  techniques  operate  to 
perform  the  same  function.   Since  analytical  comparison  is 
not  feasible,  this  functional,  comparative  approach  is 
thought  to  be  the  best  means  of  gaining  insight  into  this 
interesting  problem. 

The  purpose  of  the  analysis  is  twofold.   First,  it  is 
intended  to  indicate  those  factors  which  have  significant 
effects  on  the  performance  of  each  algorithm.   The  factors 
to  be  considered  are  environmental  variables  and  processing 
parameters.   These  are  discussed  in  Chapter  III.   Secondly, 
the  analysis  is  intended  to  point  out  the  relative  strengths 
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and  weaknesses  of  the  methods  by  comparison  of  results  on 
similar  data. 

Each  algorithm  is  evaluated  for  a  range  of  simulated 
processing  conditions.   Quantitative  and  qualitative  criteria 
are  specified  which  provide  a  comprehensive  description  of 
the  manner  in  which  each  signal  is  affected  by  processing 
for  multiple  removal.   These  criteria  also  serve,  as  a  basis 
for  comparison  of  results.   The  scope  of  the  analysis  and 
the  specific  performance  criteria  are  discussed  thoroughly 
in  Chapter  III. 
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CHAPTER  II 

ANALYTICAL  FORMULATION  AND  IMPLEMENTATION 

A *   Analytical  Formulation  of  the  TDL  Filter 

A  simple  feedback  model  for  the  propagation  of  reverberat- 
ing seismic  signals  is  given  by  Baggeroer  [3]  .   Multiple 
removal  based  on  this  model  is  then  formulated  as  an  inverse 
filtering  problem.   The  dereverberation  filter  is  designed 
using  a  least  squared  error  criterion  and  the  constraint  that 
the  filter  have  a  tapped  delay  line  structure. 

The  formulation  will  be  developed  here  from  a  different 
point  of  view  using  Baggeroer 's  feedback  model  as  a  starting 
point.   A  summary  of  the  feedback  model  is  included  for  clarity. 

Figure  3  shows  the  Laplace  transform  representation  of  a 
propagating  seismic  signal.   S(s)  is  the  transform  of  the 
source  signature.   The  signal  first  encounters  the  downward 
travel  time  delay  which  corresponds  to  a  phase  shift  in  this 
domain.   Hv.(s)  represents  the  transfer  function  of  the  earth 
beneath  the  water  column  including  the  reflections  from  layer 
boundaries  which  comprise  the  desired  information.   Internal 
multiples  or  reverberations  between  the  various  earth  layers 
are  also  included  in  H,(s).   Another  phase  shift  corresponds 
the  return  of  the  reflected  signal  through  the  water  column. 
P(s)  is  the  feedback  gain  representing  the  water  column  multiple 
mechanism.   In  most  cases  this  is  well  approximated  by  -1, 
which  corresponds  to  the  nearly  perfect  pressure  release 
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ref lection  at  the  surface.   H  (s)  includes  the  observation 

effects  such  as  hydrophone  bandwidth  and  array  tow  depth. 

Ambient  noise  and  reciever  front  end  noise  are  modelled  as 

additive  white  Gaussian  noise. 

The  overall  transfer  function  is 

-2sT 
R  (s)    H  (s)  H,  (s)  e    W 
_° =  _2 b 

S(s)     1  -  P(s)  Hb(s)  e  z    1w 

where  T   is  the  one-way  water  travel  time  and  R  (s)  is  the 

w  2  o 

received  signal  without  additive  noise. 

It  is  apparent  that  the  presence  of  multiples  is  due  only 
to  the  denominator  of  this  expression.   Thusfar  we  have 
assumed  implicitly  that  the  earth  response  can  be  modelled 
accurately  as  a  linear  system  and  that  the  multiples  are 
exactly  periodic.   The  validity  of  these  assumptions  will 
become  apparent  in  the  discussion  of  performance  in  Chapter  IV. 

The  obvious  task  is  now  to  design  an  inverse  filter  having 

the  form 

-2sT 
F(s)  =  [1  -  P(s)  Hb(s)  e    W]  . 

Hence,  we  are  required  to  estimate  T   and  the  impulse  response, 

W 

h,  (t),  corresponding  to  H,  (s).  The  earth  response  need  not 
be  estimated  precisely  for  its  entire  duration.  Estimating 
the  dominant  energy  part  of  h,(t)  is  adequate  to  produce  an 
effective  dereverberation  filter.   A  typical  deep  water 
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seismogram  has  the  great  majority  of  its  energy  concentrated 
in  the  first  200-300  msec  of  its  duration.   Effective  represent- 
ation of  this  portion  of  the  signal  requires  about  10-20 
filter  coefficients,  depending  on  the  bottom  and  source 
characteristics . 

The  transfer  function  of  equation  (1)  can  be  re-written 
in  series  form  as 


R    (s)  °°  ._  -2nsT 

O  TT  /  \  V  /  T     \    n+l         TT  /  \  W 

=    H     (s)     I        (-1)  HvJs)     e     • 

S(s)  °         n=l 

The  received  signal  then  has  the  time  domain  representation 


n+l 


r  (t)  =  s(t)  *  [h  ft)  *  I    (-1)"  x  h,  (t  -  2nT  )] 


n=l 


w 


where  *  represents  convolution.   This  can  be  rewritten  as  the 
sum  of  the  primary  return  and  the  multiple  signal. 


rQ(t)  = 


s(t)  *  hQ(t) 


n+l. 


h,  (t  -  2T  )  +  T  (-1)"  ±hK(t-2nT  ) 
b        w     L  _    '     b       w 
n=2 


rQ(t)  =  b(t)  +  m(t) 


where 


b(t)  =  s(t)  *  h  (t)  *  h^Ct-2T  ) 

o        b      w 


is  the  received  primary  and 


m(t)  =  s(t)  *  h    (t)  *  I      h,  (t-2nT) 


n=2 


w 
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is  the  received  multiple  signal. 

The  estimation  problem,  given  the  feedback  model  of 
figure  3,  is  to  determine  b(t)  in  the  presence  of  m(t)  and 
white  noise,  V7(t).   It  is  convenient  to  group  the  unwanted 
signal  components. 


n(t)  =  m(t)  +  w(t) 


(2) 


Since  the  unwanted  component  is  an  additive  one,  we 
can  consider  estimating  n(t)  and  subtracting  the  result  from 
the  received  signal.  We    then  have  the  filtering  problem 
depicted  in  figure  4 . 


r(t)  =  b(t)  +  n(t) 


f  (t) 


n  (t) 


V 


■>(+)-->  b(t) 


Figure  4 


Here  f(t)  is  the  filter  impulse  response  and  n(t)  is  the 
minimum  mean  square  error  (MMSE)  estimate  of  n(t),  given  the 
received  signal  r(t)  . 

The  number  of  digital  filter  coefficients  to  be  estimated 
is  2TW+1,  where  T  is  the  effective  duration  (portion  contain- 
ing about  80%  or  more  of  the  signal  energy)  of  h,  (t)  and  W 
is  the  signal  bandwidth.   The  coefficients  will  then  be 
spaced  at  the  Nyquist  sampling  interval  of  1/2W  seconds. 
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The  optimum  digital  filter  for  n(t)  will  have  coefficients, 
f,  ,  which  minimize 


r~ 


e  =  E  {  I        n(i/2W)  -  n(i/2W) 


1=1 


n(i/2W)  I2 


J 


(3a) 


o 


v/aere 


n(i/2W)  =  I        f  -rfii^i 

k=i    K   ^  2W 
o 


-  2T 


w 


(3D) 


The  input  in  (3b)  is  shifted  by  the  two-way  travel  time  to  avoid 
useless  filtering  of  the  signal  prior  to  the  bottom  reflection, 
[i  /2W,  if/2W]  is  the  time  interval  over  which  n(t)  is  observed. 
Substituting  (3b)  into  (3a)  yields 


e   =   E    \     I 

ti»iL 

*>   o 


-  "f 


rr  c      f(    (i-k)    n  ) 
I    .  fk  f   ~2W    '  2Tw 


k=i 


(i/2W) 


Minimizing, 


3e 

3f 


=  E 


f 
-  i  =  i 


2  I    V  r 


f   (i-k) 


2W     2Tw 


-n(i/2W) 


(i-k)    _m 
-r    OT7    -  2T 
2W        w 


0  =  5;   f.  R    (i-k)/2W 


k=i 


'k   rr 


-  R   (k/2W+2T  ) . 

nr        w 


(4) 


Here  we  have  assumed  stationarity  over  the  duration  of  the 
multiple  period.   This  assumption  has  led  to  effective  process' 
ing  of  both  real  and  synthetic  data.   From  (2) 


R   (k/2W+2T  )  =  R   (k/2W+2T  )  +  R   (k/2W+2T  ) 
nr   '      w     mr   '      w     wr   '      w' 
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Here  it  is  useful  to  assume  (see  ref.£3],  p. 15)  that  for  shifts 

close  to  2T   the  cross-correlation  of  w(t)  and  r(t)  is  very 
w 

small  compared  to  R   (t)  which  will  generally  have  a  peak  in 

this  region.   This  is  equivalent  to  assuming  that  the  white 

noise  has  a  very  short  correlation  time  compared  to  2T  .   We  then 

w 

have 


R   (k/2W+2T  )ft?  R   (k/2W+2T  ) 
nr        w     mr        w 


so  that  (4)  becomes 
i 


f 

J   f,  R   f(i-k)/2w)  =  R   (k/2W+2T  )  (4a) 

L  .   k   rr  *■      '  J  mr        w 


k=i 
o 


Baggeroer  has  derived  equation  (4a)  by  designing  the 

Wiener  filter  for  b(t)  with  the  constraint  that  the  filter  have 

a  tapped  delay  line  structure,  i.e. 

2TW 
f(t)  =  6(t)  -  I      f,  6(t-k/2W-2T  )  (5) 

k=o   k 

When  our  estimate,  n(t) ,  is  subtracted  from  the  unshifted, 
received  signal,  the  resulting  overall  filter  operation  has 
exactly  the  form  of  (5) .   This  indirect  approach  yields  the 
estimator  equations  of  reference  [3]'  without  imposing  the  TDL 
structure  directly. 

The  above  derivation  also  emphasizes  the  estimator- 
subtractor  or  prediction  error  structure  of  this  filter.   The 
entire  impulse  response  may  be  written  as  follows :  • 
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2T   zeros 
w 


This  is  the  prediction  error  structure  for  a  prediction  distance 

of  2T  .   Equation  (4a) ,  however,  which  generates  the  {f.} 

differs  from  the  prediction  equations  in  that  the  right  hand 

side  vector  is  R   (t+2T  )  rather  than  R   (t+2T  ) .   As  written, 

mr     w  rr     w  ' 

equation  (4a)  corresponds  to  the  Wiener  filter  which  produces 
the  MMSE  estimate  of  m(t)  with  r(t-2T  )  as  an  input.   Subtraction 
of  this  estimate  from  r(t)  whitens  only  those  spectral  compon- 
ents, which  are  due  to  the  multiple. 

It  is  simply  proved  that  the  magnitude  of  the  error  in  b(t) 
is  equal  to  that  in  the  prediction  operator. 

b(t)  =  r(t)  -  n(t) 

b(t)  =  b(t)  +  (n(t)-n(t)) 

|b(t)-b(t)  |  =  |n(t)-n(t)  | 

Thusfar  the  only  departures  from  optimum  estimation  have 

been  the  two  assumptions  of  stationarity  and  the  relative 

insignificance  of  R   (t+2T  ) .   One  further  assumption  is 

wr     w 

required  for  actual  implementation  of  the  filter.   Note  that 
R   (t+2T  )  is  a  required  input  which  is  apparently  not  measur- 
able from  the  given  data.   Baggeroer  has  observed  that  for  the 
deep  water  case,  which  is  of  primary  interest  for  this  method, 
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p   (T+2T  )<^  P.   (t  +  2T  )  . 
*Tnr     w    rr v    w 

That  is,  for  shifts  of  nearly  twice  the  water  travel  time,  the 
great  majority  of  the  crosscorrelation  energy  is  due  to  m(t). 
Hence,  the  equations  used  for  data  processing  are  given  by 
i 


f 

I    fk    Rrr((J-k)/2w)    =    Rrr(j/2W+2Tw)       j=0,l,...2TW 


k=i 

o 


Having  seen  the  analytical  formulation  of  this  inverse 
filtering  procedure  it  is  instructive  to  compare  it  with  the 
Backus  three-point  method.   Processing  actual  data  with  both 
filters  (ref.[3])  has  shown  the  Backus  filter  to  be  significantly 
less  effective.   Some  reasons  for  this  are  apparent  from  the 
foregoing  analysis. 

The  Backus  filter  is  rigidly  dependent  on  the  accuracy  of 
two  assumptions.   The  first,  the  assumption  of  strictly  periodic 
multiples,  is  violated  due  to  the  horizontal  separation  of 
source  and  receiver.   This  effect  becomes  more  severe  as  water 
depth  decreases.   Since  the  Backus  filter  is  implemented  as 
only  three,  equidistant  operator  coefficients  it  is  very  sensitive 
to  this  lack  of  periodicity.   Even  if  the  statistics  of  the 
signal  generate  very  accurate  estimates  of  the  bottom  reflec- 
tion coefficient  the  filter  structure  is  so  simple  and  rigid 
that  proper  cancellation  will  not  occur  if  the  multiples  are 
significantly  aperiodic. 
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A  second  restrictive  assumption  of  the  Backus  filter  is 
that  the  bottom  reflection  should  be  accurately  characterized 
by  a  single  reflection  coefficient.   All  of  the  statistical 
information  available  is  forced  into  a  single  parameter 
estimation  scheme.   It  is  apparent  that  such  a  filter  lacks 
flexibility  for  dealing  with  more  complicated  bottom  interaction 
mechanisms . 

The  relatively  better  performance  of  the  TDL  filter  on 
real  data  is  apparently  due  to  its  greater  inherent  flexibility. 
That  is,  the  finite  length  impulse  response,  or  prediction 
operator,  gives  the  filter  a  capability  for  removing  reverbera- 
tion effectively  in  cases  where  the  bottom  response  is  not 
accurately  modelled  as  a  weighted  impulse.   If  the  bottom  has 
a  ringing  or  smearing  effect  on  the  incident  signal  then  the 
deconvolution  operator  must  be  extended  in  time.   The  Backus 
filter,  because  of  its  rigid  structure,  cannot  accomodate  these 
situations.   The  TDL  structure  provides  2TW+1  (usually  10-20) 
parameters  v/hich  can  be  varied  in  the  design  procedure  to 
optimize  dereverberation  of  each  seismogram.   The  special  case 
of  an  ideal  bottom  will  generate  a  filter  response  which  is 
essentially  a  single  spike  proportional  to  the  bottom  reflection 
coefficient.   This  result  has  been  confirmed  in  the  analysis  of 
synthetic  data.   In  such  a  case  the  TDL  filter  consists 
basically  of  the  first  two  points  of  the  Backus  three-point 
filter.   Performance  (multiple  energy  removed)  in  these  cases 
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was  found  to  be  essentially  independent  of  operator  length. 

The  structural  flexibility  of  the  TDL  filter  also  gives 
it  some  potential  for  dealing  with  aperiodic  multiples.   It 
should  be  noted,  however,  that  the  TDL  algorithm,  like  the 
Backus  and  other  classical  dereverberation  techniques,  is 
essentially  a  correlation-cancellation  operation.   Consequently, 
increased  aperiodicity  of  multiples  can  be  expected  to  degrade 
performance  in  all  cases. 

B.   Implementation  of  the  TDL  Algorithm 

Figure  5  shows  a  flow  diagram  of  the  filter  implementation 
used  for  this  analysis.   Actual  programming  was  done  in 
Fortran  IV  for  use  on  a  32K  computer.   Referring  to  figure  5, 
the  correlation  function  is  computed  by  the  standard  shift-and- 
add  operation  with  no  windowing  applied.   Results  of  windowing 
are  included  in  reference  [3].   Correlation  time  is  variable 
and  may  be  specified  by  the  operator.   The  crosscorrelation 
function  is  approximated  by  the  correlation  function  as  discussed 
in  the  previous  section. 

Solution  of  the  filter  equations  is  accomplished  by 
conventional  matrix  inversion.   The  T°eplitz  symmetry  can  be 
exploited  for  computational  savings.   Spacing  of  the  operator 
elements  is  determined  by  the  estimated  signal  bandwidth  which 
is  specified  as  an  input  parameter. 

Actual  deconvolution  is  implemented  exactly  as  shown  in 


-24- 

RECEIVED 
SIGNAL 


CORRELATION    FCN 


R   (t) 

rr 


APPROXIMATE 
CROSSCORRELATION 
FCN 

R   (x+2T  ) 

mr      w 


SOLUTION  OF 
ESTIMATOR  EQUATIONS 

T  f .R   (T. )=R   (T.+2T  ) 
L      1  rr   1    mr   l    w 


IMPLEMENT  TDL 


b(t)=r(t)*{fi) 


DE REVERBERATED 
SIGNAL 


Figure  5   Flov;  chart  of  TDL  algorithm. 
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figure  6,  i.e.,  by  means  of  a  tapped  delay  line  or,  equivalently , 
convolution  with  the  gapped  operator. 


C.   Analytical  Formulation  of  Homomorphic  Deconvolution 

A  homomorphic  system  is  one  which  obeys  a  generalized 
principle  of  superposition  and  which  can  be  represented  as  an 
algebraically  linear  transformation  between  two  vector  spaces. 
A  detailed  description  of  the  theory  of  homomorphic  systems 
is  given  by  Oppenheim  and  Schafer  [6] .   This  material  will  not 
be  repeated  in  depth  here;  rather,  we  shall  discuss  the  basic 
characteristics  of  homomorphic  systems  for  convolution  with 
emphasis  on  those  properties  which  facilitate  dereverberation. 
Additional  discussions  of  these  properties  are  found  in 
references  [5] ,  [7]  and  [8]  . 

The  usefulness  of  linear  systems  for'  separating  additively 
combined  signals  is  due  primarily  to  superposition.   Signals 
which  are  added  and  happen  to  be  disjoint  in  the  frequency 
domain  can  be  separated  by  means  of  an  appropriate  bandpass 
filter.   Homomorphic  systems  for  convolution  have  a  similar 
effect  on  signals  which  have  been  convolved.   That  is,  a 
homomorphic  transformation  maps  the  input  signal  to  a  domain 
in  which  the  convolved  components  may  be  disjoint.   Such  a 
transformation  is  illustrated  in  figure  7. 


x(n> 


x(n) 


Figure  7 
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The  mapping  characteristic  of  this  system  is  intuitively 
apparent.   Suppose  x(n)  is  composed  of  two  components  which 
have  been  convolved. 

x(n)  =  x1 (n)  *  x2 (n) 

The  z-transform  operation  yields  a  signal  with  multiplied 
components,  X,(z)  and  X2(z).   The  logarithm  output  is 

X(z)  =  logtX^z)]  +  log  [X2(z)]. 
The  inverse  z-transform, 

/S  /N  /N 

x(n)  =  x1 (n)  +  x2 (n)    , 

preserves  the  additive  combination  of  the  components  and  yields 
a  sequence  which  is  real  and  stable  for  a  real  and  stable  x(n) . 
The  sequence  x(n)  is  called  the  complex  cepstrum  of  x(n). 
Although  it  is  real  for  real  inputs,  the  "complex"  is  retained 
to  emphasize  that  it  contains  both  the  magnitude  and  phase 
information  from  X(z).   Hence  the  complex  logarithm  is  required 
even  for  real  x(n).   (We  shall  omit  the  modifier  here  for 
brevity. ) 

The  cepstrum  variable,  n,  is  normally  called  the  quefrency 
(a  paraphrase  of  frequency)  or  period  variable.   Filter  opera- 
tions in  this  domain  are  generally  similar  to  those  encountered 
in  the  frequency  domain.   Exact  subtraction  of  x, (n)  from  x(n) 
followed  by  computation  of  the  inverse  cepstrum  (figure  8) 
yields  the  sequence  x~ (n) ,  exactly,  in  the  time  domain. 
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x(n) 


x(n) 


Figure  8 


It  is  this  property  which  renders  homomorphic  processing  valu- 
able for  deconvolution. 

As  an  example  of  complex  cepstrum  transformation  consider 
a  signal  composed  of  a  short  pulse  (2  samples)  convolved  with 
a  decaying,  periodic  impulse  train. 


1 


where 


x(n)  =  s  (n)  *  p(n) 

00 

x(n)  =  (6(n)  +  i  o(n+l))*  £  (-1)  kRkS  (n-kT) 

*  .k=o 


Rl  <  1   and  T  >  2 


Taking  the  z-transform, 


v  k  k  -kT 

X(z)  =  (1  +  z/2)  •  I     (-1TNR  z 


k=0 


-T 


-1 


X(z)  =  (1  +  z/2) • (1  +  Rz   ) 


The  logarithm  then  produces  a  sum 


X(z)  =  logCl  +  z/2)  -  logCl  +  Rz   )  . 
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Both  terms  are  simply  expanded  in  Laurent  series. 


log  U  +  z/2)  =    I 
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t-ir     -n 


n=-l  n2 
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z  <  2 


oo       n  n 

i    /t    t,  ~T    v   (-1)  R     -nT    .   -T. 
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n=l 
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The  z-transforms  are  easily  recognized. 

n 


sM  =  -i-y 


n2 


-n 


,    n  =  -1,  -2  .  .  .-<» 


r-l^kRk 
p(n)  =  l  ±j     6(n-kT)   ,  k  =  1,2,...' 


x(n)  =  s(n)  +  p(n) 


V7e  see  that  s  (n)  occupies  only  the  negative  quefrency 
region  and  p(n)  only  the  positive  quefrency  region.   Exact 
deconvolution  can  be  accomplished  in  this  case  by  zeroing  the 
desired  half  of  the  cepstrum. 

The  canonic  form  of  homomorphic  systems  for  convolution  is 
shown  in  figure  9 . 


x(n) 


y(n) 


Figure  9 
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The  characteristic  systems  D^  and  D~   are  shown  in  figures  7 
and  8  respectively.   The  system  L  is  a  conventional  linear 
system.   When  L  has  a  transfer  function  of  1,  then  x(n)  is 
recovered  exactly  at  the  output  of  D^  .   The  choice  of  L  v/ill 
determine  the  effectiveness  of  deconvolution  for  a  given  input 
sequence . 

Two  further  specifications  are  required  to  ensure  the 
validity  and  uniqueness  of  the  transformation.   The  complex 
logarithm  is  a  multivalued  function  with  an  infinite  number  of 
branches . 

log[X(z)]  =  log|x(z) |  +  jf Arg[X (z) ]±  2ukJ       for  all  k 

v/here  Arg  specifies  the  principal  value  of  the  phase.   This 
ambiguity  must  be  resolved  while  simultaneously  satisfying  the 
requirement  that  X(z)  be  a  valid  z-transform.   Note  that  if 
x(n)  is  to  be  real  and  stable,  X(z)  must  be  conjugate  symmetric 
and  analytic  in  an  annulus  of  the  z-plane  containing  the  unit 
circle.   That  is,  the  real  and  imaginary  parts  of  X(z)  must  be 
continuous  functions  of  z  in  the  region  including  the  unit 
circle.   The  imaginary  part,  arg[X(z)],  can  only  be  made 
continuous  by  "unwrapping"  Arg[x(z)]  in  such  a  way  that  all 
jumps  of  ±  2Trk  are  removed.   This  unique  unwrapping  leads  to  a 
unique,  valid  X(z)  which  transforms  to  a  stable,  real  x(n). 
The  requirement  of  a  continuous  phase  curve  poses  some  computa- 
tional difficulties  which  will  be  discussed  in  the  following 
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section . 

It  is  apparent  from  the  foregoing  description  that 
homomorphic  systems  have  the  potential  for  separating  convolved 
signals.   One  might  expect,  however,  by  analogy  with  linear 
systems  that  deconvolution  is  most  effective  for  signals  with 
certain  cepstral  properties.   This  is,  in  fact,  the  case  and, 
fortunately,  seismic  signal  components  are  generally  amenable 
to  deconvolution.   Recall  that  a  seismogram  is  modelled  as  a 
convolution  of  a  source  signature  and  an  impulsive  reflector 
series.   Reverberation  appears  as  a  minimum  phase,  periodic 
addition  to  the  reflector  series.   Thus,  we  have 

rQ(n)  =  p(n)  *  (b(n)  +  m(n)] 

where  p(n)  is  the  source  signature,  b(n)  is  the  desired  signal 
and  m(n)  is  reverberation.   Generalizations  can  be  made  concern- 
ing the  cepstral  properties  of  each  component. 

The  source  signature  is,  in  general,  a  mixed  phase,  short 
duration  time  sequence.   It  is  clear  from  the  definition  of  the 
z-transform  that  any  such  finite  sequence  transforms  to  a 
rational  function  of  z  with  no  poles  except  at  the  origin.   In 
general 

,  m        -,  I 
P(z)  =  C  z~K  n  Cl-a.z"x)  n  (1-b.z)    |a.|,|b  |<  1 
i=l     x     j=l     J  J 

The  a.  and  b.  represent  zeros  inside  and  outside  the  unit  circle 
i       D 

respectively.   z~   corresponds  to  a  linear  phase  shift. 
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m  .    % 
P(z)  =  log  C  +  I    log(l-a.z   )  +  £  log(l-b.z) 
i-1      1      j=l      J 


log  C         n=0 
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m 

a. 

I 
i=l 
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I 

b.n 

I 
j  =  l 

n 

n>°      .  (6) 


n<0 


Here  it  has  been  assumed  that  the  linear  phase  term  is  removed 
before  computation.   p(n)  is  a  two-sided  sequence  which  is 
always  of  infinite  duration  but  decays  faster  than  1/n . 
Hence,  most  of  the  cepstral  energy  is  concentrated  near  the 
quefrency  origin. 

The  reflector  sequence  is  modelled  as  a  train  of  randomly 
spaced  impulses  which  may  be  mixed  phase.   Stoffa,  Buhl  and 
Bryan  [5]  give  a  general,  but  very  complicated  expression 
for  the  complex  cepstrum  of  such  a  sequence.   Some  specific 
examples  are  given  by  Schafer  [7].   The  resulting  cepstrum  is 
an  impulse  train  with  impulses  at  the  time  domain  impulse 
locations,  at  all  their  multiples,  and  at  various  other  loca- 
tions, both  positive  and  negative  on  the  quefrency  axis.   Three 
important  observations  can  be  made. 

The  cepstrum  of  a  minimum  phase  reflector  train  contains 
no  contributions  for  negative  quefrency.   Consider  the  special 
case  of  equation  (6)-  in  which  all  the  b.  are  zero.   This 
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corresponds  to  all  zeros  being  inside  the  unit  circle,  i.e.,  a 
minimum  phase  sequence.   Note  that  p(n)  and  s (n)  in  the  example 
are  minimum  and  maximum  phase  respectively,  leading  to  causal 
and  anticausal  cepstra. 

Secondly,  it  happens  that  no  non-zero  cepstral  contributions 
occur  between  the  origin  (first  impulse)  and  the  location  of 
the  second  impulse  in  time,  if  the  time  series  is  minimum  phase. 
Therefore,  the  cepstrum  of  a  minimum  phase  impulse  train,  unlike 
that  of  a  general  sequence  does  not  have  its  contributions 
concentrated  near  the  origin.   The  cepstrum  of  a  minimum  phase 
impulse  train  will  always  contain  a  gap  equal  to  that  between 
the  first  two  time  domain  contributions. 

Finally,  we  observe  (see  Schafer  [7])  that  a  reflector 
series  can  easily  be  made  minimum  phase  by  exponential  weighting. 

r"  (n)  =  wnr(n)         |w|  <  1 

R'  (z)  =  C  z"k  n  fl-(aiw)z~1j  fl-(bj/w)zj 

The  value  of  w  is  chosen  so  that 

I b .w   I  >  1      for  all  b . . 

1  J    '  D 

Weighting  of  the  impulse  train  is  effected  by  weighting  of  the 
entire  signal  since,  if 

s  (n)  =  p(n)  *  b(n)   , 
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then 

wns(n)  =  (wnp(n))*(wnb(n)). 

Note  that  the  reflector  series  may  be  made  minimum  phase  without 
making  the  entire  signal  minimum  phase.   Very  little  weighting 
is  normally  required  in  practice. 

The  remaining  contribution  to  the  seismic  signal  is  reverbera- 
tion.  This  component  is  merely  a  special  case  of  a  minimum 
phase  impulse  train  in  which  the  impulses  are  periodic.   It  is 
easily  verified  that  if 

a 

m(n)  =  I    y  (n)6  (n-kn  ) 
k=l  ° 

then 

m(n)  =  y  (n/nQ)  . 

The  cepstrum  is,  therefore,  periodic  with  the  same  period  as 
the  reverberation. 

A  useful  property  for  deverberation  is  derived  by  Stoffa, 
et  al.[5].   The  derivation  is  summarized  here  because  of  its 
direct  pertinence  to  seismic  processing. 

Consider  a  normalized  multiple  signal, 

00 

m(n)  =  I    (-l)1R16(n-i2T  )  (7) 

i=0 

where  R  is  the  bottom  reflection  coefficient.   Then,  as  we 

have  seen 

i^i 


m(n)  =  I       C""11  R   S(n-2iT  ). 
i=l 


!  W 
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Subtracting  the  first j  cepstral  contributions  and  transforming, 

m.  (n)  =  mCn)  -  7   C"1)  R   6  (n-2iT  ) 
3  L  .•  w' 

i=l    1 

°°     Di   ~2iT 

M.tz)  =  I        -  £-2- 

3      .  .      i 
i=j+l    x 

OO 

M. (z)  =  n     exp  [--  z"2lTwJ  . 
3  i=j+l       i 

Expanding  in  a  power  series 

00    «>   _ik  -2ikT 

m.(z)  -  n   I  -   5-z w 

J     i=j+l  k=0   kli 

00  ,  ^  -2T  ,i+k  oo   oo        -2T   2D+^ 

m.(z)  =i+l  il55 !£  +  j  £   J^ s> 

D         k-l     3+k     k=3  i=1  2<i+j)(j+k-i) 

(8) 

Comparing  (7)  and  (8),  the  first  j  time  domain  multiples  have 
been  removed  completely  and  the  (j+l)st  multiple  is  reduced  by 
l/(j+l).    All  succeeding  multiples  are  also  reduced.   Thus, 
removal   of  only  the  first  cepstrum  multiple  would  remove  the 
first  time  domain  multiple  and  reduce  the  second  by  1/2,  the 
third  by  1/3,  etc.. 

Having  seen  the  cepstral  properties  of  each  seismic  signal 
component  the  advantages  of  homomorphic  deconvolution  are 
apparent.   The  source  signature  and  reflector  series  have  their 
cepstral  energy  concentrated  in  different  regions  of  the  quefrency 
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axis,  thus  facilitating  removal  of  the  source.   The  cepstral 
contributions  of  the  reverberation  occur  at  predictable  locat- 
tions  so  that  multiple  removal  is  possible. 

Due  to  the  nature  of  the  homomorphic  transformation, 
linear  filtering  of  the  cepstrum  is  not  the  normal  convolution 
operation  but  a  simple  zeroing  of  the  unwanted  contributions. 
That  is,  the  linear  filter  is  generally  frequency  invariant 
rather  than  time  or  quefrency  invariant.   The  name  "quefrency" 
was  adopted  to  reflect  this  reversal  of  the  customary  time  and 
frequency  filtering  roles. 

The  foregoing  analysis  is  based  on  assumptions  similar  to 
those  employed  in  classical  seismic  processing.   Namely, 
we  have  assumed  that  the  seismogram  consists  of  a  source 
signature  convolved  with  distinct,  impulsive  reflectors  and 
periodic  multiples.   It  is  difficult  to  predict  the  sensitivity 
of  the  overall  processing  scheme  to  these  assumptions  because 
of  its  complex  analytical  structure.   Hence,  various  parameters 
have  been  varied  in  the  performance  analysis  to  obtain  an 
empirical  measure  of  this  sensitivity. 

Finally,  we  note  that  the  additive  noise  was  not  included 
in  the  analytical  formulation  of  the  homomorphic  processing 
scheme.   The  algorithm  is  designed  to  separate  convolved 
components  and,  unfortunately,  no  effective  processing  gain  is 
achieved  over  added  noise.   In  practive,  additive  noise  has 
been  dealt  with  through  classical  bandpass  filtering.   This 


-37- 

performance  analysis  includes  a  description  of  the  effects  of 
additive  noise  on  homomorphic  deconvolution. 

D.   Implementation  of  Homomorphic  Deconvolution 

Figures  7 ,  8  and  9  show  the  sequence  of  operations  required 
for  homomorphic  deconvolution.   A  processing  scheme  was  designed 
to  implement  this  algorithm  in  Fortran  IV  on  a  32K  digital 
computer.   A  flow  diagram  of  the  scheme  is  shown  in  figure  10. 
Some  aspects  of  the  computation  are  noteworthy. 

All  z-transforms  in  the  algorithm  are  implemented  via  FFT. 
Recall  from  the  analytical  formulations  that  z-transforms 
involved  in  the  processing  of  a  real,  stable  sequence  are 
required  to  have  regions  of  convergence  which  include  the  unit 
circle.   The  discrete  Fourier  transform  is  simply  a  sampling 
of  the  z-transform  on  the  unit  circle  which,  for  properly  band- 
limited  signals,  is  sufficient  to  specify  the  signal  completely. 
Data  sequences  are  normally  padded  with  zeros  to  reduce  cepstral 
aliasing,  e.g.,  2048-point  cepstra  are  computed  for  1024-point 
seismograms.   Since  the  cepstrum  is  always  of  infinite  extent, 
a  truncated  version  always  results  in  some  aliasing  when 
computation  is  not  done  recursively. 

The  major  difficulty  in  computing  an  accurate  cepstrum  is 
the  computation  of  a  continuous  phase  curve.   The    ta  sequences 
are  normally  sampled  at  a  rate  based  on  the  frequency  content. 
There  is  no  assurance,  however,  that  this  sampling  rate  is  adequate 
to  uniquely  specify  x(z)  =  log[X(z)]. 
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The  method  used  for  this  analysis  is  due  to  Tribolet  [9]  and 
is  thought  to  be  quite  accurate  and  efficient.   The  flow 

diagram  for  this  algorithm  is  shown  in  figure  11. 

dXX  z) 
The  phase  derivative,  — ■—— —  and  principal  value,  Arg[X(z)], 

are  easily  computed  from  the  transforms  of  x(n)  and  n  x(n) 

(see  appendix  A) .   The  values  of  the  derivative  are  integrated 

using  the  trapezoidal  rule  and  the  integral  output  is  compared 

with  Arg  X(z)  at  each  step.   If  the  two  values  do  not  agree 

within 

2niT  ±  e        n  =  0,±1,  ±2,  .  .  . 

where  e  is  a  small  positive  number,  the  latest  computed  value 
of  arg  X(z) ,  say  a. ,    is  discarded. -  The  routine  then  returns 
to  the  last  correct  integration  value,  computes  an  intermediate 
derivative  value,  and  begins  integrating  with  a  step  size  half 
that  of  the  original  grid.   The  integrate-and-compare  process 
is  continued  at  this  step  size  until  a.  has  been  computed 
correctly  or  until  a  comparison  fails.   Integration  is  resumed 
at  the  initial  step  size  in  the  former  case  or,  in  the  latter, 
step  size  is  again  halved.   The  number  of  possible  step  sizes 
is  theoretically  unlimited.   The  value  of  e   may  be  adjusted 
by  trial-and-error  for  most  efficient  integration.   The 
adaptive  step  feature  compensates  for  the  under sampling  problem 
in  a  very  efficient  and  accurate  manner. 

Linear  phase  contributions  are  easily  identified  and 
removed  from  the  computed  continuous  phase  curve.   Having 
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computed  the  log  magnitude  and  phase,  the  inverse  z-transform 
is  straightforward. 

Linear  filtering  is  accomplished  by  zeroing  the  unwanted 
cepstrum  values.   One  might  consider  windowing  procedures  which 
are  common  in  linear  filtering,  but  these  were  not  employed  in 
the  present  analysis. 

The  inverse  cepstrum  computation  is  completely  straight- 
forward since  no  ambiguities  arise  in  the  exponentiation 
process.   The  final  steps  are  shifting  the  output  sequence 
by  the  linear  phase  value  and  unweighting  the  shifted  sequence, 
if  necessary. 
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CHAPTER  III 
DESCRIPTION  OF  THE  PERFORMANCE  ANALYSIS 

The  purpose  of  this  analysis  is  to  evaluate,  both  quantita- 
tively and  qualitatively,  the  performance  of  the  TDL  and  homo- 
mo  rphic  dereverberation  techniques.   Each  method  is  to  be 
evaluated  for  a  variety  of  simulated  seismic  conditions  in  order 
to  determine  those  factors  which  significantly  influence  perfor- 
mance.  The  comparative  nature  of  the  analysis  is  intended  to 
emphasize  the  relative  strengths  and  weaknesses  of  each  technique 

It  should  be  noted  that  absolute  performance  figures  are 
not  inherently  valuable,  especially  when  obtained  from  synthetic 
data.   The  diverse  geological  and  oceanographic  conditions 
encountered  in  marine  seismology  coupled  with  the  many  different 
processing  systems  currently  employed  may  be  expected  to  yield 
a  range  of  absolute  results.   The  greater  value  of  this  analysis 
is  to  indicate  the  parameters,  environmental  and  mathematical, 
which  can  be  expected  to  affect  significantly  the  performance 
of  these  algorithms.   The  numbers  obtained  provide  a  measure 
of  the  relative  performance  of  the  two  methods  under  similar 
conditions  and,  in  some  cases,  provide  asymptotic  performance 
criteria.   Synthetic  data  were  chosen  so  that  signal  parameters 
could  be  accurately  controlled. 

Three  criteria  are  specified  for  comprehensive  evaluation 
of  performance.   The  first,  most  direct,  measure  of  effectiveness 
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is  the  percentage  of  multiple  energy  removed  from  the  signal. 
This  is  easily  measured  by  calculating  the  zero-lag  correlation 
of  the  signal,  in  time  windows  spanning  only  the  multiple 
locations,  before  and  after  processing.   That  is,  the  squared 
amplitude  (energy)  of  the  signal  in  the  multiple  region  is 
computed  after  processing  and  subtracted  from  the  squared 
amplitude  of  the  original  signal  in  that  region.   This 
difference  is  divided  by  the  original  energy  of  the  multiple 
to  yield  the  fraction  of  energy  removed  by  processing.   The 
use  of  synthetic  data  facilitates  this  method  of  analysis 
since  reflectors  and  multiples  can  be  placed  in  disjoint  regions 
to  avoid  ambiguity. 

The  second  criterion  is  the  amount  of  signal  distortion 
introduced  by  dereverberation.   This  is  measured  by  comparing 
reflector  energy  before  and  after  processing.   As  before, 
reflectors  and  multiples  must  be  disjoint  for  meaningful  results. 

Although  these  two  criteria  provide  an  accurate  quantitative 
measure  of  performance  they  are  restricted  to  situations  in 
which  reflectors  and  multiples  do  not  overlap.   The  overlap 
case  is  most  important  in  processing  real  data  since  the 
multiples  then  abscure  the  reflectors  most  severely.   In  order 
to  judge  performance  in  these  situations  we  must  evaluate 
qualitatively  the  improvement  of  visual  interpretability . 
This  visual  enhancement  of  reflector-to-multiple  ratio  is  our 
third  criterion.   It  is  an  important  measure  in  spite  of  its 
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subjective  nature  since  the  primary  means  of  seismogram 
analysis  is  visual  interpretation. 

The  analysis  is  limited  in  scope  to  the  single  channel 
processing  configuration.   Although  both  techniques  are 
applicable,  in  principle,  to  multichannel  processing,  the  many 
additional  variables  involved  and  unavailability  of  appropriate 
synthetic  data  would  lead  to  a  complicated  extension  of  this 
analysis.   A  single  channel  treatment  is  adequate  to  identify 
the  important  performance  traits  of  both  methods. 

Parameters  to  be  varied  fall  into  two  general  categories; 
environmental  and  operational.   The  environomental  parameters 
include  noise  level,  multiple  periodicity,  multiple-to-signal 
level  and  multiple/reflector  separation  or  overlap.   These  are 
varied  within  ranges  which  are  thought  -to  be  representative  of 
ambient   conditions  normally  encountered  in  marine  seismology. 
Effects  of  noise  level  are  considered  only  for  the  case  of 
white  Gaussian  noise.   The  effects  of  aperiodicity  have  not 
been  evaluated  for  the  TDL  algorithm  because  it  is  intended 
primarily  for  deep  water  use  where  only  the  first  multiple  is 
usually  of  interest.   Operational  parameters  refer   to  those 
which  can  be  controlled  during  processing.   These  .include  filter 
cutoff  frequencies,  operator  lengths,  travel  time  estimates, 
cepstral  stopbands,  and  weighting.   Windowing  of  the  correlation 
function  is  not  evaluated  although  a  discussion  of  this  subject 
is  contained  in  reference  [3] . 
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The  algorithm  used  in  generating  simulated  earth  impulse 
responses  is  due  to  Theriault  [10J .   A  brief  description  of 
Theriault's  earth  model  is  given  here. 

The  model  is  based  on  the  following  assumptions: 
CI)  An  air  gun  source  generates  longitudinal  pressure 
waves  which  impinge  upon  all  earth  layers  as 
normally  incident  plane  waves. 

(2)  All  earth  layers  are  horizontally  infinite, 
parallel  and  homogeneous. 

(3)  Abrupt  changes  in  acoustic  impedance  occur  at 
each  layer  interface  and  these  boundaries  are 
characterized  by  the  customary  acoustic  reflec- 
ion  and  transmission  coefficients. 

(4)  The  earth  has  a  uniform  density. 

(5)  The  water  column  is  a  non-attenuating  fluid 
with  a  perfect  pressure  release  interface  at 
its  surface. 

(6)  Each  layer  is  characterized  by  a  transfer 
function,  F(w)  which  represents  the  attenuation 
and  travel  time  delay  for  that  layer. 

These  assumptions  are  incorporated  into  a  lumped  parameter 
model.   Figure  12  shows  a  frequency  domain  model  of  a  two 
layer  earth.   The  F.  (co)  have  the  functional  form 
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Figure  12   Tv;o  layer  earth  transfer  function  schematic 
(after  Theriault) . 
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where  xi ,  c.  and  ou  are  the  ith  layer  thickness,  sound  speed 
and  layer  attenuation  parameter  respectively.   These  transfer 
functions  are  combined  using  a  semi-group  property  rather  than 
the  usual  frequency  domain  multiplication. 

The  -1  multiplier  completes  a  feedback  loop  around  the 
source   which  generates  the  water  column  multiples.   Since  the 
water  column  is  assumed  to  have  no  attenuation  the  multiples 
appear  in  the  earth  response  as  impulses,  and  in  the  resulting 
seismogram  as  replicas  of  the  source  signature  reduced  by  the 
bottom  reflection  coefficient. 

The  above  multiple  mechanism  is  inadequate  for  representing 
effects  of  incoherent  (distorted)  multiples  and  varying  bottom 
interaction  mechanisms.   These  effects  are  introduced  by  inser- 
tion of  water  column  attenuation  which  causes  the  bottom  response 
and  multiple  responses  to  be  of  exponential  form  given  by  the 
Fourier  transform  of  (9) .   Thus  the  bottom  response  is  extended 
in  time  and   each  multiple  is  a  distorted  version  of  the  previous 
one.   Examples  of  multiple  appearance  with  and  without  water 
column  attenuation  are  shown  in  figure  13a  and  b. 

The  topmost  loop  of  figure  12  simulates  the  effects  of 
finite  receiver  depth  t.   Any  number  of  layers  with  the  desired 
parameters  can  be  combined  into  an  overall  earth  transfer 
function  which  is  easily  transformed  to  yield  the  earth  impulse 
response . 

Synthetic  seismograras  for  this  analysis  were  generated  by 
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Figure  13   (a)  Synthetic  sei  sinogram  with  coherent  multiples, 
(b)  Synthetic  sei'^onr^  with  distorted  multiples 
due  to  water  column  attenuation. 
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convolution  of  various  Theriault  impulse  responses  with  air 
gun  signatures  obtained  from  actual  at-sea  recordings.   A 
typical  signature  is  shov/n  in  figure  14. 

It  is  important  to  note  that  these  seismograms  contain 
all  water  column  multiples  and  internal  multiples  of  all 
orders.   For  realistic  parameter  values  earth  attenuation 
characteristics  usually  render  internal  multiples  negligible. 

Finally,  we  note  one  drawback  of  using  this  model  for 
multiple  removal  analysis.   The  algorithm  produces  multiples 
which  are  exactly  periodic.   This  periodicity  gives  these 
seismograms  strong  correlation  characteristics  which  are  not 
usually  encountered  in  practice.   The  lack  of  periodicity  in 
actual  seismograms  is  due  to  the  horizontal  separation  of  the 
source  and  receiving  array.   The  difference  in  travel  paths 
arising  from  this  separation  is  illustrated  below  for  a  primary 
return  and  first  multiple. 
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Figure  14  Air  gun  signature  (obtained,  by  at-sea  recordings) 
convolved  with  synthetic  earth  response  functions 
to  form  seismograms . 
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For  a  known  separation  and  water  depth  the  travel  time 
difference  can  be  easily  calculated.   This  effect  becomes 
minimal  in  deep  water  where  cos  a  and  cos  3  are  approximately 
equal  to  one . 
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CHAPTER  IV 
RESULTS 

A.   Introduction 

The  results  of  applying  the  multiple  removal  methods 
described  in  Chapter  II  to  synthetic  data  are  described  here  in 
terms  of  the  criteria  of  Chapter  III.   Performance  based  on  the 
first  two  criteria  is  emphasized  because  it  can  be  quantified 
much  more  accurately.   Specifically,  the  reflectors  and 
multiples  have  been  positioned  at  distinct  locations  in  most 
cases  so  that  the  amount  of  energy  removed  from  reflectors  and 
multiples  can  be  measured  without  ambiguity.   This  emphasis 
leads,  however,  to  a  certain  lack  of  realism  in  several  of  the 
synthetic  data  plots.   Separation  of  this  kind  in  an  actual 
seismogram  would,  of  course,  eliminate  the  necessity  for 
multiple  removal.   Hence,  several  cases  of  interfering  reflectors 
and  multiples  are  also  shown.   Although  these  are  not  amenable 
to  quantitative  analysis  they  can  be  judged  on  the  basis  of 
the  third  criterion,  viz.,  improvement  of  visual  record  quality. 

A  second  deviation  from  normal  processing  conditions  has 
been  required  to  compare  effectively  the  performance  of  the 
two  methods.   Homomorphic  dereverber'ation  is  accomplished  in 
practice  (see  {5J )  concurrently  with  source  deconvolution, 
where  practical,  since  both  operations  are  simply  performed 
after  the  cepstrum  has  been  computed.   This  leads  to  a 
considerable  amount  of  energy  removal  at  each  multiple  and 
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ref lector  location  which  is  due  to  source  deconvolution  alone. 
Although  this  may  lead  to  improvement  of  the  record,  the  criterion 
of  energy  removed  does  not  accurately  measure  dereverberation 
performance  in  the  same  way  it  does  for  the  TDL  results.   For 
example,  source  deconvolution  might  typically  lead  to  90% 
reduction  of  multiples  and  90%  removal  of  reflector  energy. 
This  gross  change  in  configuration  cannot  be  effectively 
compared  with  TDL  processing  on  the  basis  of  energy  removed. 
Hence,  source  removal  has  not  been  accomplished  in  most  cases 
tested.   The  results  of  this  "partial"  processing  can  be 
quantitatively  evaluated  and  easily  compared  with  TDL  results. 
Cepstral  filtering  which  includes  source  removal,  thus  retain- 
ing only  high  quefrency  energy,  is  referred  to  as  "longpass 
filtering".   Several  examples  of  longpass  filtering  are  included 
and  interpreted  in  terms  of  the  third  criterion. 

The  performance  of  each  method  is  discussed  individually 
for  various  processing  situations.   A  direct  comparison  of  the 
performance  of  both  methods  for  the  same  data  is  also  included. 
The  comparison  is  extended  in  Chapter  V. 

B.   Results  of  TDL  Dereverberation  Performance 

1.  Operator  Length,  Multiple  Distortion  and  Multiple-to- 
Signal  Ratio 

The  number  of  tap  gains  required  for  optimum  multiple 

removal  was  found  to  be  highly  signal-dependent.   Recall  from 

Chapter  II  that  the  tapped  delay  line  is  essentially  an 
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estimate  of  the  high  energy  portion  of  h,(t),  the  earth  impulse 
response.   In  most  applications  the  reflected  bottom  response 
is  dominant.   The  time-bandwidth  product  of  this  response 
would  then  be  expected  to  govern  the  filter  length  requirements. 
Measured  results  confirm  this. 

Seismograms  containing  exactly  impulsive  (one  sample) 
bottom  responses  exhibited  little  or  no  variation  of  performance 
with  filter  length.   Typical  filter  impulse  responses  for  a 
seismogram  of  this  type  are  shown  in  figure  15.   The  reflection 
coefficient  in  this  case  is  0.3.   The  first  tap  gain  is  close  to 
-0.3  in  each  response,  as  would  be  expected  from  the  Backus 
formulation  in  which  the  second  operator  point  is  an  estimate 
of  the  bottom  reflection  coefficient.   Increasing  filter  length 
can  be  seen  to  cause  variation  in  the  "estimation"  of  the  reflec- 
tion coefficient.   Figure  16  shows  energy  removed  vs.  filter 
length  for  this  seismogram.   Multiple  energy  removed  decreases 
slightly  with  increasing  filter  length.   Reflector  distortion 
is  nearly  constant  at  low  values  (6%  for  one  and  -2%  for  the 
other) .   The  signal  used  in  figure  16  is  shown  in  figure  17a. 
The  processed  result  shown  was  obtained  vising  only  one  tap. 

Introduction  of  a  non-impulsive  multiple  mechanism  was 
found  to  produce  a  marked  dependence  of  performance  on  operator 
length.   Synthetic  seismograms  were  generated  with  finite 
length  bottom  responses,  resulting  in  distorted  multiples 
which  are  not  simply  weighted  replicas  of  the  source  signature. 
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Figure  15   TDL  impulse  responses  for  three  operator  lengths 
(a)  5  taps  (b)  10  taps  (c)  15  taps 
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Figure  16   Operator  length  vs.  performance  for  a  signal  with 
with  impulsive  multiples. 
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Figure  17   (a)  Seismogram  with  coherent  multiples. 

(b)  Result  of  TDL  processing  v/ith  a  one  point 
operator. 
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Figure  18  shows  filter  performance  vs.  operator  length 
for  three  seismograms  having  different  degrees  of  distortion 
in  their  bottom  reflection  mechanisms.   This  distortion  is 
equivalent  to  extension  of  the  bottom  impulse  response  in  time. 
The  signals  associated  with  curves  (1) ,  12)    and  C3)  are   shown 
in  figures  19a,  19b  and  19c  respectively.   The  increasing 
distortion  of  the  multiple  at  2.0  seconds  is  evident.   Curve  (1) 
corresponds  to  a  signal  with  very  slight  multiple  distortion. 
Operator  length   is  seen  to  have  no  effect  on  performance. 
The  seismogram  corresponding  to  curve  (2)  has  a  bottom  impulse 
response  which  is  significant  for  T  =  .03  seconds.   The  tap 
spacing  in  this  case  is  1/2W  =  4.8  8  msec,  so  that  2TW  =  6.1, 
and  six  or  seven  tap  gains  should  be  adequate  if  the  signal 
has  been  properly  sampled.   Reference  to  curve  (2)  confirms 
that  increasing  the  filter  length  beyond  seven  does  not  improve 
performance.   Filters  of  fewer  than  seven  elements  yield 
monotonically  decreasing  performance.    The  bottom  response 
for  curve  (3)  is  significant  for  .045  seconds  so  that,  for  the 
same  tap  spacing,  2TW  =  9.2  and  we  anticipate  that  nine  or  ten 
taps  will  be  adequate.   This  is,  in  fact,,  the  case. 

Figure  20  shows  the  5,  9  and  15  point  filter  impulse 
responses  associated  with  curve  (3) .   There  is  an  observable 

0     —  l-..f 

convergence  to  a  t  e~    shape  which  is  the  actual  functional 
form  of  the  synthetic  bottom  response.   Figures  20a  and  20b 
exhibit  a  "diverging  tail"  effect  which  was  found  to  be  common 
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Figure  18   Performance  vs.  operator  length  for  three  signals 
with  bottom  interaction  times  varying  from  (1) 
impulsive  to  (3)  .045  seconds. 
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Figure  19   Seismograns  associated  with  the  curves  of  figure 

18.   (a)  Curve  (1).   (b)  Curve  (2).   (c)  Curve  (3). 
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Figure  2  0   Filter  impulse  responses  associated  with  Curve  (3) 
of  figure  18.   (a)  5  points   (b)  9  points 
(c)  15  points. 
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when  the  specified  filter  length  is  too  short.   This  type  of 
phenomenon  occurs  in  some  numerical  approximation  methods  when 
an  inadequate  number  of  terms  is  specified. 

It  is  apparent  from  the  above  results  that  the  effect  of 
operator  length  is  closely  related  to  the  bottom  interaction 
mechanism.   As  discussed  in  Chapter  II,  the  filter  should  be 
an  estimated  replica  of  the  bottom  impulse  response  when  the 
interaction  process  can  be  accurately  modelled  as  a  convolution 
of  the  source  signature  with  the  bottom  response.   Figures 
15  and  20  are  good  examples  of  this  behavior. 

In  the  cases  summarized  in  figure  18  performance  increases 
as  multiple  distortion  increases.   This  need  not  be  true  in 
general  since  bottom  interactions  may  become  very  complex. 
The  slowly  varying  t  e     responses  lead  to  operators  which 
have  a  greater  cancellation  effect  as  they  are  extended.   A 
higher  bandwidth  bottom  impulse  response  might  not  exhibit 
this  behavior.  As  it  was  not  possible  to  include  more  complicated 
bottom  responses  in  the  earth  model  used,  these  effects  were 
not  investigated  further. 

Reflector  distortion  was  found  to  be  nearly  constant  for 
all  filter  lengths  tested.   The  reflector  at  2 . 7  seconds  was 
essentially  undistorted  in  all  cases.   The  3.5  second  reflector 
had  an  average  distortion  of  7%.   Figures  21-24  show  some 
processed  results  for  the  signals  in  figure  19.   Each  of  the 
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Figure  21   (a)  Seismogram  with  very  coherent  multiples. 
(b)  Result  of  TDL  processing  with  three  taps 
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Figure  22   (a)  Seismogram  with  moderately  distorted  multiples; 

bottom  response  is  significant  for  .035  seconds, 
(b)  Result  of  processing  with  7  taps. 
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Figure  2  3   (a)  Seis^ogram  with  considerable  multiple  distortion; 

bottom  response  is  significant  for  .045  seconds, 
(b)  Result  of  TDL  processing  with  11  taps. 
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first  three  figures  (21,  22  and  23)  show  signals  before  and 
after  processing  with  the  fewest  number  of  taps  required  to 
achieve  optimum  performance.   Visually,  multiple  removal  is 
almost  complete  in  each  case.   The  signal  of  figure  23a  is 
shown  in  figure  2  4  after  processing  with  3,  5  and  7  taps. 
The  improvement  from  figure  24a  to  24c  is  obvious. 

Multiple-to-signal  ratio  (defined  here  as  the  ratio  of 
energy  in  the  first  multiple  to  that  in  the  largest  sub-bottom 
reflector,  abbreviated  MSR)  was  found  to  be  of  little  importance 
in  most  cases  of  interest.   For  signals  in  which  multiples  are 
large  enough  to  be  a  problem  (comparable  to,  or  larger  than 
smaller  reflectors) ,  the  multiple  dominate  the  crosscorrelation 
function  so  that  an  effective  filter  is  generated.   For  these 
cases  the  performance  was  found  to  be  insensitive  to  the  width 
of  the  time  window  used  for  correlation.   In  the  relatively 
less  interesting  case  of  signals  with  small  multiples,  the 
performance  is  greatly  dependent  on  the  choice  of  correlation 
window.   If  large  reflectors  are  included  in  the  window,  perfor- 
mance is  adversely  affected  because  of  the  large  reflector 
contribution  to  the  statistics.   For  some  cases  of  interest 
this  effect  may  be  a  consideration  in  choosing  an  appropriate 
correlation  window.   Inclusion  of  large  reflectors  should  be 
avoided. 
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Figure  2  4   Results  of  processing  seismogram  of  figure  23a  with 

inadequate  TDL  lengths  (a)  3  taps  (b)  5  taps  (c)  7  taps 
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2.  Water  Travel  Time  Estimate 

It  was  noted  in  Chapter  U  .that  an  estimate  of  the  water 

column  travel  time  is  required  for  implementation  of  the  TDL 

algorithm.   This  is  accomplished  in  practice  by  various  methods 

including  visual  estimation,  energy  detection  and  correlation 

techniques.   The  estimate  appears  in  the  filter  design  equations 

as  the  minimum  shift  of  the  crosscorrelation  function.  R 

mr 

This  estimate  may  also  be  identified  as  the  prediction  distance 
when  the  operator  is  interpreted  as  a  prediction  error  filter. 

The  actual  estimation  of  water  column  travel  time  was  not 
investigated  in  this  analysis.   The  effects  of  travel  time 
estimation  on  filter  performance  were,  however,  considered. 

Figure  25  shows  the  results  of  water  travel  time  estimation 
errors  for  the  ideal  case  of  a  signal  with  impulsive  multiples. 
The  unprocessed  seismogram,  shown  in  figure  26a,  contains 
reflectors  at  2.7  and  3.5  seconds  and  a  strong  multiple  at  2.0 
seconds.   The  actual  two-way  travel  time  in  this  case  is  1.0 
second.   The  strong  similarity  between  the  bottom  reflection 
and  multiple  is  apparent. 

First  multiple  energy  removed  is  very  sensitive  to  travel 
time   estimation,  with  an  error  of  ±  10  msec  resulting  in  a 
performance  degradation  of  about  20%.   This  effect  is  analogous 
to  that  observed  in  matched  filter  receivers  in  that  the 
coherent  signal  components  exhibit  very  high  correlation  for  a 
small  range  of  lags.   In  this  case  the  strong  coherence  and 
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Figure  25   Effects  of  v;ater  column  travel  time  estimate  on 

multiple  removed  for  a  signal  with  coherent  multiples. 
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Figure    2  6       (a)    Unprocessed   seisraogram  with   impulsive   multiples. 
(b)    TDL   result   based   on    correct  water   travel    time 
estimation . 
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Figure  26  (cont'd)  (c)  Result  of  TDL  processing  based  on  an 

early  water  travel  tine  estimate. 
(c)  Same  processing  based  on  a  late  estimate. 
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multiple  periodicity  yield  the  sharp  correlation  peak  for  lags 
near  the  water  column  travel  time,  which  is  also  the  multiple 
period. 

The  second  and  third  multiple  performance  peaks  occur 
at  1.08  and  1.04  seconds  respectively,  and  have  values  of  24% 
and  60%.   These  shifted  and  reduced  peaks  are  caused  by 
digitization  effects.   For  such  an  extremely  coherent  signal, 
a  multiple  position  deviation  of  one  sample  (due  to  sampling 
interval  round-off)  can  cause  the  cancellation  operation 
to  be  severely  affected.   In  this  case  the  second  and  third 
multiple  performance  maxima  are  due  to  secondary  crosscorrelation 
peaks  introduced  by  the  periodic  oscillations  in  the  source 
signature.   These  peaks  occur  at  intervals  approximately 
equal  to  the  period  of  the  basic  source  (air  gun)  frequency. 
The  three  performance  peak  values  correspond  closely  to  the  first 
three  air  gun  pulse  amplitudes. 

Figures  26  b,  c  and  d  show  the  results  of  on-time,  early 
and  late  travel  time  estimates  respectively.   In  the  first  case 
the  multiple  has  been  effectively  removed  while  the  early  and 
late  estimates  lead  to  multiples  which  are  still  significant 
after  processing. 

The  effects  of  travel  time  estimation  error  on  the 
operator  impulse  response  are  seen  in  figure  27.   The  optimum 
estimate  yields  an  operator  with  a  large  peak  at  the  origin, 
nearly  equal  to  the  bottom  reflection  coefficient   (.2),  and  a 
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Figure  2  7   TDL  filter  impulse  responses  based  on  travel  time 

estimates  which  are  (a)  40  msec  early,  (b)  correct, 
and  (c)  20  msec  late. 
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smaller  peak  at  the  source  pulse  oscillation  interval.   The 
initial  spike  is  due  to  the  largest  peak  in  the  crosscorrelation 
function  Capproximated  here  by  the  correlation  function)  which 
occurs  at  a  shift  equal  to  the  two-way  travel  time.   The  smaller 
peak  in  the  filter  is  introduced  by  an  additional  shift  of  one 
source  period.   The  early  (by  40  msec)  estimate  still  allows 
observation  of  the  crosscorrelation  maximum  and  yields  the 
shifted  spike  of  figure  27a.   It  is  evident  in  figure  25  that 
performance  is  reasonably  good  for  estimation  errors  less  than 
about  .07  seconds,  which  is  the  length  of  the  tapped  delay 
line.   For  greater  travel  time  errors  the  range  of  correlation 
lags'  computed  does  not  include  the  peak;  hence,  the  filter  is 
ineffective.   Late  travel  time  always  results  in  poor  perfor- 
mance since  the  peak  is  not  observed.   Figure  2  7c  shows  a 
typical  impulse  response  due  to  late  travel  time  estimation. 
The  seismogram  of  figure  2  8a  contains  a  small  reflector 
at  2.2  seconds  and  a  larger  one  at  3 . 0  seconds.   The  first 
multiple  partially  overlaps  the  smaller  reflector  and  the  second 
multiple  coincides  with  the  larger.   Figures  28  b,  c  and  d 
demonstrate  how  travel  time   estimation  errors  can  lead  to 
ambiguous  results.   Figure  28b  was  processed  using  the  correct 
travel  time,  resulting  in  good  resolution  of  both  reflectors. 
The  early  and  late  estimates  lead  to  the  signals  of  figures 
28c  and  d  in  which  the  smaller  reflector  cannot  be  clearly 
resolved. 
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Figure  28   (a)  Seismogram  with  overlapping  reflectors  and 

multiples . 
(b)  Result  of  TDL  processing  based  on   a  correct 
travel  time  estimate. 
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Figure    2  8    (cont'd) 


(c)  Result  of  TDL  processing  based  on  a 
40  msec  early  travel  time  estimate. 

(d)  Same  processing  based  on  a  20  msec  late 
estimate . 
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Less  coherent  seismograms  exhibit  less  sensitivity  to 
travel  time  estimation  as  shown  in  figure  29.   Curves  CD  , 
(2)  and  (3)  were  generated  using  the  seismograms  of  figure  19  a, 
b  and  c  respectively,  which  contain  increasing  distortion  in 
their  multiples.   The  MSR  is  considerably  lower  than  in 
figure  25  so  that  the  later  multiples  are  not  clearly  visible. 
The  crosscorrelation  peaks  are  considerably  broader  than  in 
the  impulsive  multiple  signals.   Curve  (3),  the  least  coherent 
signal,  exhibits  the  lowest  sensitivity  to  travel  time  estimates. 
Increasingly  sharper  peaks  are  evident  for  curves  (2)  and  (1) . 
The  performance  degradation  due  to  over-estimation  is  precipitous 
in  all  cases  although  the  rate  of  fall-off  is  related  to 
signal  coherence. 

Reflector  distortion  (not  shown)  was  negligible  in  all 
cases  for  the  first  reflector  and  nearly  constant  at  7%  for 
the  second. 

Examples  of  processed  signals  for  curves  (2)  and  (3) 
are  shown  in  figures  3  0  and  31. 

3.  Multiple  Periodicity 

The  effects  of  aperiodic  multiples  on  TDL  dereverberation 
performance  were  not  investigated  in  this  analysis  since  the 
algorithm  is  designed  primarily  for  deep  water  use,  where  only 
one  multiple  is  significant  in  most  applications.   Successful 
employment  of  this  technique  in  shallow  water  is  limited  since 
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Figure  29   Effects  of  water  column  travel  time  on  three  signals 
with  varying  degrees  of  multiple  distortion. 


-79- 


Rl 


R2 


' 


b; 


B2 


B3 


I 


A/w— 


JU/_ 


fVW~ Am- 


(a) 


Vw~- 


(b) 


1.0 


2.0 


3.0 


4.0 


5.0 


Figure  30   Processed  seismograms  associated  with  figure  19b 
and  curve  (2)  of  figure  29.   (a)  Travel  time 
estimate  40  msec  early  (b)  40  msec  late. 
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Figure  31   Processed  seismograns  associated  with  figure  19c  and 
curve  (3)  of  figure  29.   (a)  Travel  tine  estimate  40 
.msec  early  (b)  40  msec  late. 
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the  approximation  of  the  crosscorrelation,  R   ,  by  the 
correlation,  R   ,  near  the  multiple  onset  is  not  generally 
valid.   Reflector  energy  in  the  multiple  regions  is  usually 
significant  in  shallow  water  seismograms . 

It  is  apparent  from  the  TDL  formulation  and  the  discussions 
of  the  preceding   two  sections  that  the  TDL  algorithm  has 
some  capability  to  compensate  for  aperiodicity  in  which  the 
deviation  is  of  the  order  of  the  tapped  delay  line  length. 
For  greater  aperiodicities  the  TDL  filter  will  not  be  effective 
on  later  multiples.   Figure  25  indicates  that  even  slight 
deviations  can  be  very  detrimental  in  later  multiple  removal. 
Dynamic  corrections  can  be  applied  as  suggested  by  Backus  [1] , 
however,  these  are  beyond  the  scope  of  this  study. 

4.   Reflector/Multiple  Overlap 

Figures  32-36  show  seismograms  before  and  after  processing 
for  several  cases  in  which  multiples  and  reflectors  overlap. 
A  brief  description  of  each  situation  is  given  here. 

The  seismogram  of  figure  32a  has  reflectors  at  1.5,  2.1 
and  3.0  seconds,  and  distorted  multiples.   Significant  energy 
near  2.0  seconds  is  removed  by  processing  but  a  clearly  visible 
response  is  still  present.   Other  regions  of  the  signal  are 
not  visibly  affected.   The  visual  resolution  is  not  significantly 
improved  in  this  case,  however,  stacking  of  successive  shots 
after  multiple  removal  could  be  employed  to  enhance  the  reflector 
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Figure  32   Sei sinogram  with  reflectors  at  1.5,  2.1  and  3.0 

seconds  (a)  before  and  (b)  after  TDL  processing 
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in  this  situation.   The  success  of  stacking  would  require 
that  the  reflector  be  more  coherent  than  the  residual  multiple 
after  dereverberation. 

In  figure  33  reflectors  occur  at  1.5,  2.02  and  3.0  seconds. 
Reflectors  may  be  expected  to  have  a  strong  influence  on  the 
crosscorrelation  function  in  this  situation  of  extreme  overlap. 
The  result  is  nearly  complete  removal  of  the  second  reflector 
with  the  multiple.   The  last  reflector  is  clearly  visible  at 
3.0  seconds.   In  this  example  the  tapped  delay  line  length 
(.054  second)  is  greater  than  the  reflector/multiple  separation 
(.02  second) .   The  bottom  impulse  response  is  significant  for 
.035  second  so  that  it  cannot  be  estimated  accurately  with  an 
operator  which  is  shorter  than  the  reflector/multiple  separation. 
Even  if  the  bottom  response  were  much  shorter,  the  effect  of 
the  large  reflector  in  the  crosscorrelation  window  would  lead 
to  degraded  performance.   This  example  shows  worse  degradation 
than  would  be  expected  in  practice  because  the  reflectors  are 
very  coherent  in  this  synthetic  data.   The  less  coherent 
reflections  in  deep  water  signals  would  generally  be  less 
susceptable  to  removal.   When  reflector  energy  is  significant 
in  the  crosscorrelation,  however,  and  TDL  length  is  greater  than 
reflector/multiple  separation,  the  filter  design  algorithm 
essentially  interprets  the  reflector  as  part  of  the  multiple 
and  tries  to  remove  it. 

Figure  34  contains  reflectors  at  1.6,  1.95  and  3.0  seconds. 
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Figure  33   Seismogram  with  reflectors  at  1.5,  2.02  and  3.0 
seconds  (a)  before  and  (b)  after  TDL  processing 
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Figure  34   Seismogram  with,  reflectors  at  1.6,  1.95  and  3.0 
seconds  (a)  before  and  (b)  after  TDL  processing 
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These  locations  are  slightly  different  from  those  of  figure  33. 
Significant  energy  is  present  at  the  second  reflector  location 
after  processing.   Since  the  reflector  precedes  the  first 
multiple  in  this  case  it  has  a  less  significant  effect  on  the 
crosscorrelation  function  than  the  2.02  second  reflector  in  the 
previous  example.   The  variability  of  performance  on  individual 
seismograms  (of  similar  structure)  with  overlap  again  suggests 
the  use  of  stacking  to  enhance  reflectors. 

The  reflectors  of  figure  35  are  situated  as  in  figure  34; 
however,  the  first  reflector  and  multiple  are  significantly  larger 
in  this  case.   This  seismogram  is  more  typical  of  a  shallow 
water  return.   The  processing  is  quite  effective  on  the  first 
and  second  multiples.   The  aperiodicity  of  reverberation  in 
most  shallow  water  seismograms  would  lead  to  much  poorer 
performance  on  later  multiples,  in  practice. 

Figure  36  illustrates  a  situation  where  the  overlap  is 
not  severe  but  visual  resolution  is  impaired  by  the  first  and 
second  multiples.   Reflectors  occur  at  1.55,  2.4  and  3.25 
seconds.   Processing  results  in  effective  reduction  of  both 
multiples  and  significantly  improved  resolution  in  the  second 
and  third  reflectors. 

We  conclude  from  these  examples  of  ref lector/multiple 
overlap  that  visual  improvement  due  to  dereverberation  varies 
widely,  depending-  upon  several  aspects  of  the  individual  signal 
structure.   The  bottom  impulse  response,  the  extent  of  reflector/ 


-87- 


I 


Rl    R2 


B2,  R3 


Bl 


/y\A^ — ^/^j/i/wVW^^ — ^-av 


(a) 


B3 


wWv^- 


(b) 


1.0 


2.0 


3.0        4.0 


5.0 


Figure  35   Seismograra  with  reflectors  at  1.6,  1.95  and  3.0 
seconds  (a)  before  and  (b)  after  TDL  processing 
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Figure  36   Seismogram  with  reflectors  at  1.55,  2.4  and  3.25 
seconds  (a)  before  and  (b)  after  TDL  processing 
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multiple  separation,  the  relative  coherence  of  reflector  and 
multiples,  and  the  relative  sizes  of  signal  components  all 
appear  to  affect  performance.   Each  of  these  factors  influences 
the  effectiveness  of  the  crosscorrelation  operation  in  identify- 
ing energy  which  is  specifically  due  to  the  multiples. 

The  visual  improvement  due  to  dereverberation  can  best  be 
determined  by  analysis  of  continuous  seismic  profiles  (presenta- 
tions showing  many  seismograms  side-by-side)  since  coherent 
residual  energy,  if  present,  becomes  apparent  as  visual  trends 
in  the  data.   The  effects  of  stacking  can  also  be  observed  in 
this  format.   It  was  not  practical  to  generate  such  a  profile 
with  synthetic  data  but  the  single-shot  results  indicate  that 
TDL  dereverberation  is  potentially  useful  for  enhancing  the 
visibility  of  reflectors  which  are  partially  masked  by 
multiples . 

5.   Additive  White  Noise  and  Filtering 

White  Gaussian  noise  was  generated  in  the  following  manner. 
First,  a  set  of  uniformly  distributed  random  numbers  was 
obtained  using  a  standard  digital  routine.   A  set  of  approximately 
Gaussian  numbers  was  then  formed  by  summing  separate  groups  of 
twelve  of  the  uniformly  distributed  numbers.   The  resulting  set 
was  weighted  to  obtain  the  desired  standard  deviation.   Figure 
37  shows  a  seismogram  with  two  different  levels  of  added  noise, 
each  lowpass  filtered  at  50  Hz. 
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Figure  37   Seismograms  with  added  white  noise,  lowpass  at  50  Hz. 
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Figure  38  illustrates  the  performance  of  the  TDL  filter 
for  various  noise  levels  on  a  seismogram  having  disjoint 
reflectors  and  multiples.   The  maximum  noise  standard  deviation 
shown  corresponds  to  a  SNR  of  0.7  dB.   No  bandpass  filtering 
was  done  prior  to  dereverberation.   The  results  presented  here 
are  typical  of  the  cases  tested.   Multiple  energy  removed 
decreases  monotonically  with  increasing  noise  level;  almost 
linearly  in  the  case  of  the  first  multiple.   The  later  multiples 
are  more  severely  affected.   As  the  noise  level  becomes  comparable 
to  their  amplitudes,  the  filter  becomes  ineffective  in  removing 
them.   This  behavior  is  due  to  the  effect  of  increasing  incoherent 
energy  in  the  signal.   The  correlation-cancellation  operation 
is  designed  to  detect  coherent,  periodic  components.   These 
become  increasingly  masked  as  more  noise  is  added.   For  this 
reason  seismograms  with  larger  multiples  exhibit  less  sensitivity 
to  noise.   Seismograms  of  similar  structure  (i.e.,  exactly  the 
same  reflector  and  multiple  locations)  whose  multiple  energies 
differed  by  a  factor  of  twenty  were  found  to  exhibit  consider- 
able dereverberation  performance  differences  in  high  noise. 
Processing  of  the  signal  with  the  larger  .multiples  resulted  in 
25%  greater  removal  of  first  multiple  energy. 

Lowpass  filtering  before  dereverberation  can  lead  to  a 
significant  improvement  in  TDL  filter  performance  in  some 
signals.   The  plotted  results  of  lowpass  filtering  two  noisy 
signals  at  various  frequencies  are  shown  in  figure  39. 
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Figure  3  8   Effects  of  noise  level  on  TDL  performance  for  a 
signal  with  coherent  periodic  multiples. 
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Figure  39   TDL  performance  vs.  lowpass  filter  cutoff  frequency 

for  signals  with  coherent  (lower  curve)  and  distorted 
(upper  curve)  multiples. 
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CA  third  order  Butterworth  filter  was  used  throughout.   The 
visual   difference  in  two  filtered  signals  is  shown  in  figure  40. 
The  first  is  lowpass  filtered  at  90  Hz  and  the  second  at  30  Hz.) 
The  seismogram  used  to  generate  the  upper  curve  contains 
multiples  which  are  considerably  distorted  while  the  lower 
curve  is  due  to  a  more  coherent  signal.   Decreasing  the  pass- 
band  from  90  to  30  Hz  produces  a  12%  improvement  in  performance 
in  the  second  case  but  the  less  coherent  signal  is  relatively 
insensitive  to  filtering.   Several  other   signals  evaluated 
were  found  to  exhibit  the  same  behavior.   The  phenomenon  is 
apparently  due  to  the  averaging,  or  smoothing  effect  of  the 
operator  in  the  distortion  case.   Recall  that  signals  of  this 
kind  tend  to  have  more  extended  waveforms  in  their  TDL  operators 
(See  figure  20).  Convolution  of  such  functions  with  a  noisy 
signal  "smoothes  out"  the  noise  because  of  the  significant 
operator  length  and  the  random  fluctuation  of  the  noise.   This 
has  the  twofold  result  of  better  overall  performance  in  noise 
and  less  sensitivity  to  removal  of  higher  frequency  noise  energy. 

Performance  degenerates  for  filter  cutoff  frequencies  below 
30  Hz  because  the  spectral  energy  of  the. signal  itself  is 
concentrated  in  this  range.   Figure  41  shows  the  distribution 
of  energy  in  the  frequency  domain  for  the  seismogram  associated 
with  the  upper  curve  of  figure  39.   Most  of  the  energy  is 
concentrated  between  5  and  30  Hz. 

Figures  42-45  illustrate  the  visual  improvement  of  noisy, 
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Figure  40   Seismograms  with  equal  white  noise  levels,  lowpass 
filtered  at  (a)  90  Hz  and  (b)  30  Hz. 
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Figure  41   Squared  magnitude  of  the  Fourier  transform  for  the 
signal  associated  with  the  upper  curve  in  figure  39 
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Figure  42   (a)  Unprocessed  seismogram  with  added  noise,  Fc=50  Hz. 
(b)  Result  of  TDL  processing. 
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Figure  43  (a)  Seismogram  of  figure  42a  with  high  noise,  Fc=50  Hz- 
(b)  Result  of  TDL  processing. 
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Figure  44   (a)  Seismogram  with  low  noise,  F  =50  Hz. 
(b)  Result  of  TDL  processing. 
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Figure  45   (a)  Seismogram  of  figure  44  with  high  noise,  F  =50  Hz. 
(b)  Result  of  TDL  processing. 
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lowpass  filtered  C50Hz)  signals  after  processing.   Two  different 
seisraograms  are  shown,  each  at  two  different  noise  levels, 
before  and  after  processing.   It  can  be  seen  from  these  figures 
that  no  serious  reflector  distortion  occurs  due  to  noise  or 
filtering. 

The  presence  of  white  noise  has  a  predictable  effect  on 
the  appearance  of  the  TDL  operator.   The  impulse  response  itself 
becomes  noisy  due  to  the  added  incoherent  energy  and  a  bias  is 
introduced  into  the  large  cancellation  peak.   An  example  of  this 
is  shown  in  figure  46  for  a  signal  with  coherent  multiples. 
The  variation  in  the  first  tap  gain  is  significant  and  the  change 
in  appearance  of  the  overall  impulse  response  is  appearent. 

C.   Results  of  Homomorphic  Dereverberation  Performance 

1.  Introduction 

The  theory  of  homomorphic  dereverberation  is  based  on  the 
properties  of  periodic  minimum  phase  impulse  train  cepstra.   The 
incorporation  of  more  realistic  effects,  such  as  aperiodicity , 
distorted  multiples,  and  non-impulsive  reflectors,  leads  to 
analytical  intractability  in  most  cases.   Hence,  in  presenting 
the  experimental  results,  we  view  the  more  complex  signal 
structures  in  terms  of  their  relation  to  the  ideal  signal  models 
of  Chapter  II.   In  so  doing  we  hope  to  provide  a  basic  reference 
for  interpreting  observed  phenomena,  and  to  provide  a  stronger 
intuitive  picture  of  the  mechanisms  at  work  in  homomorphic 
dereverberation . 
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All  cepstra  used  in  this  analysis  are  2048  points  in 
length.   Seismograms  contain  1024  points  except  in  cases  where 
resampling  was  required.   Zeros  have  been  appended  as  necessary. 
The  2048  point  length  has  been  used  throughout  the  analysis  for 
uniformity  and  reduction  of  aliasing.   In  most  practical  seismic 
processing,  shorter  cepstrum  lengths  are  adequate. 

2.   Multiple  Periodicity 

It  was  shown  in  Chapter  II  that  later  components  of  a 
periodic,  minimum  phase  impulse  train  can  be  significantly 
reduced  by  removal  of  the  first  one  or  two  cepstral  contributions 
This  property  is  potentially  valuable  for  dereverberation , 
especially  in  shallow  water  cases  where  several  strong  multiples 
may  appear  in  the  data.   In  order  for  this  property  to  be  useful 
it  must  be  relatively  insensitive  to  (at  least)  slight 
aperiodicity  since  actual  reverberation  is  not  exactly  periodic. 
The  result  of  S toff a,  et  al.  summarized   in  Chapter  II  is 
extended  here  to  rapidly  decaying,  aperiodic,  minimum  phase 
impulse  trains. 

Consider  removing  the  first  j  cepstral  contributions  of 
mCn)  ,  the  cepstrum  of  a  minimum  phase  impulse  train  with 
contributions  C-R)   at  n.,  i  =  1,2,3,...,  and  |r|  <  1. 


We  then  have 


m  In)  =  mCn)  -  \   i-|i-  6  (n-n^l 
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Trans forming , 
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Exponentiating , 
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This  may  be   expressed   as    a  power   series. 
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The  rapidly  decaying  coefficients  of  z    &  combine  multiplica- 
tively  to  yield  the  time  domain  impulse  coefficients  of  m.(n). 
If  only  the  first  cepstral  impulse  is  removed  (j=l) ,  the 
largest  value  attained  by  the  third  coefficient  in  brackets  is 


C-R) 
1? 


21 
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which  is  insignificant  for  reflection  coefficients  of  interest. 
All  succeeding  terms  in  CIO)  will  be  vanishingly  small.   We 
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then  have  the  approximate  expression 


m.  Cz)  =  n 

3  A=j+1 


Expanding  this  product  we  obtain 


lJ 


M  (z)  =  1  +  J   ££L     -  z  ^k  +  l  R)       z   23+k+ 

3  k=l  -1™  k=3  K 


which  transforms  to 

f-R}j  +  1  (--r)3+2 

m.  (n)  =  5(n)  +  -i— ^-  .   6(11-11..,)  +  -^-^-  «  6  (n-n  ._)  +  ..  . 
3  j+1  J+l     j  +  2  3+2 

(ID 
It  is  apparent  that  the  remaining  terms  have  been  reduced  by 

factors  equal  to  those  obtained  in  the  periodic  case. 

The  above  result  has  been  verified  experimentally  by  process- 
ing a  periodic  signal  containing  only  a  strong  bottom  reflection 
(R=.55)  and  multiples  (figure  47a).   Removal  of  the  first 
multiple  component  in  the  cepstrum  results  in  75%  and  85%  energy 
reductions  of  the  second  and  third  multiples  respectively. 
These  correspond  exactly  to  the  50%  and  67%  amplitude  reductions 
predicted  by  (11) .   The  periodicity  was  then  disturbed  by 
average  deviations  of  1%,  5%,  10%  and  20%  of  the  original 
period.   In  each  case  the  later  multiple  reductions  coincided 
with  (11) .   Figure  47  shows  waveforms  before  and  after  process- 
ing for  the  exactly  periodic  signal  and  the  case  of  20% 
aperiodicity . 
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Figure  47   (a)  Signal  composed  of  only  periodic  multiples  of 

the  bottom  response. 
(b)  Result  of  removing  the  first  multiple  cepstrum 
contribution. 
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Figure  47  (cont'd) 


(c)  Signal  of  (a)  with  periodicity  disturbed 
by  average  deviations  of  20%  (.2  sec). 

(d)  Result  of  removing  the  first  multiple 
cepstrum  contribution. 


-108- 

3.   Multiple-to-Signal  Ratio 

The  results  obtained  for  aperiodic  reverberation  must  be 
extended  to  signals  containing  reflectors  as  well  as  multiples 
if  they  are  to  be  of  practical  use.   This  is  not  feasible 
analytically  due  to  the  nonlinearity  of  the  homomorphic  trans- 
formation and  the  reflector  characteristics  in  the  earth  model 
used.   Recall  that,  for  this  analysis,  reflectors  have  the 
form  t  e    rather  than  simply  impulses.   Each  seismogram  con- 
tains an  additive  combination  of  these  reflectors  with 
multiples.   The  logarithmic  operation  on  the  Fourier  transform 
of  this  sum  causes  the  cepstral  contributions  due  to  reflectors 
and  multiples  to  be  analytically  indistinguishable.   Hence, 
this  phenomenon  has  been  investigated  empirically  for  various 
reflector  sizes.   It  was  found  that  cepstral  properties  of 
multiples  are  essentially  preserved  in  the  presence  of  non- 
impulsive  reflectors  although  important  effects  were  evident 
in  the  cases  tested. 

Percentage  removal  of  the  first  multiple  was  found  to  be 
dependent  on  the  width  of  the  cepstral  stopband.   Figure  4  8 
Illustrates  this  effect  for  seismograms  with  different  reflector 
stengths  and  periodic  multiples.   Each  signal  requires  a 
stopband  of  about  200  msec  (41  samples)  for  complete  removal 
of  the  first  multiple.   As  notch  width  is  decreased  the  perfor- 
mance becomes  increasingly  sensitive  to  MSR.   For  the  smallest 
notch  width  shown  (40  msec) ,  performance  varies  almost  30% 
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with  MSR.   In  the  absence  of  reflectors  the  first  multiple 
could  be  completely  removed  by  zeroing  one  sample  of  the  cepstrum, 
It  appears  that  the  inclusion  of  reflectors  has  the  effect  of 
distributing  the  multiple  energy  in  the  quefrency  domain, 
although,  in  all  cases  tested  the  energy  remains  concentrated 
near  the  time  domain  multiple  location.   Figure  49  shows  the 
0.5  second  section  of  the  cepstrum  centered  at  the  multiple 
location  for  each  of  the  seismograms  of  figure  48.   The  cepstra 
are  identical  in  form  but  the  absolute  cepstral  energy  increases 
with  MSR.   Three  large  peaks  occur  in  each  cepstrum  between 
.94  and  1.0  second.   This  similarity  in  form  suggests  that 
equal  stopbands  should  produce  equivalent  percentage  reduction 
of  multiples.   We  conclude,  however,  from  the  experimental 
results  that  a  greater  portion  of  the  multiple  energy  is 
concentrated  near  1.0  second  in  the  higher  MSR  cepstra.   This 
effect  makes  it  difficult  to  select  stopband  limits  by  peak 
detection  or  visual  inspection  of  the  cepstrum.   Notch  widths 
which  yield  the  best  trade-off  between  multiple  reduction  and 
reflector  distortion  must  be  determined  by  trial-and-error 
for  particular  applications. 

4.  Multiple  Distortion 

Since  homomorphic  dereverberation  is  theoretically  based 
on  a  model  of  strictly  impulsive  multiples,  it  is  important  to 
observe  the  performance  of  this  technique  in  the  more  realistic 
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case  of  an  extended  bottom  interaction  mechanism.   We  employ 
the  same  bottom  response  functional  form  discussed  earlier 
in  this  chapter.   The  seismogram  of  figure  50a  contains 
distorted  multiples  due  to  a  bottom  response  of  .045  seconds 
duration.   Reflectors  are  present  at  2 . 7  and  3.5  seconds. 
Figure  50b  shows  the  region  of  the  cepstrum  centered  at  the 
first  multiple  location.   The  peaks  are  much  broader  than  those 
of  figure  49.   The  dominant  energy  is  concentrated  near  the 
multiple  location  as  before  but  it  now  appears  to  be  more 
distributed. 

Table  1  summarizes  the  effects  of  varying  the  stopband 
width  and  location  in  the  cepstrum  of  figure  50b.   First  and 
third  multiple  energy  removed  varies  widely  while  the  second 
multiple  is  not  significantly  affected. by  the  passband  dimensions 
In  some  cases,  extending  the  stopband  decreases  the  amount  of 
multiple  energy  removed.   Zeroing  the  .94-. 9 8  second  region 
results  in  greatly  increased  reduction  of  the  third  multiple 
(74-78%)  but  first  multiple  reduction  is  degraded  by  10-12%. 
This  behavior,  which  has  been  observed  in  several  cases,  is 
not  predicted  by  the  theory  and  is  thought  to  be  a  computational 
effect.   The  high  amplitude  oscillations  which  are  dominant  in 
the  left  side  C. 75-1.0  second)  of  figure  50b,  or  the  low 
frequency  "drift"  which  is  apparent  throughout  the  section  may 
be  computational  noise  which  contributes  to  this  phenomenon. 

Several  observations  can  be  made  in  this  case.   Although 
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Figure  50   (a)  Seismogram  with  distorted  multiples 
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Stopband   Limits  ENERGY   REMOVED 

( S6C ) 

1st  MULT   2nd  MULT   3rd  MULT   1st  REFL.  2nd  REFL. 


.94-1.06 

.75 

.48 

.78 

-.13 

-.01 

.96-1.06 

.78 

.44 

.74 

-.13 

.01 

.96-1.00 

.65 

.49 

.56 

-.12 

-.03 

.98-1.02 

.79 

.44 

.45 

-.12 

-.03 

.98-1.04 

.88 

.44 

.43 

-.12 

-.03 

.98-1.06 

.89 

.45 

.54 

-.» 

-.02 

.99-1.04 

.78 

.42 

.39 

-.12 

-.05 

1.0-1.04 

.64 

.45 

.48 

-.13 

-.06 

1.02-1.06 

.50 

.45 

.50 

-.13 

1 
-.06 

i 

TABLE  1 


EFFECTS  OF  STOPBAND  LIMITS  ON  PERFORMANCE 
FOR  THE  SIGNAL  OF  FIGURE  50. 
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performance  is  somewhat  erratic,  a  stopband  of  40-80  msec 
which  straddles  the  first  multiple  onset  time  yields  significant 
multiple  reduction.   The  peak  performance  in  this  case  is  about 
10%  below  that  achieved  with  undistorted  multiples  but  the 
sensitivity  to  minor  stopband  variations  is  less  dramatic. 
Reflector  distortion  does  not  increase  significantly  due  to  the 
presence  of  distorted  multiples. 

Figure  51  shows  processed  seismograms  resulting  from  the 
application  of  various  stopbands  to  the  cepstrum  of  figure  50. 

Finally,  we  note  that  the  effects  of  aperiodicity  have 
not  been  discussed  in  the  distorted  multiple  case.   Due  to 
limitations  of  the  earth  model  it  was  not  possible  to  investigate 
these  effects.   Thusfar,  however,  we  have  relaxed  the  theoretical 
assumptions  regarding  periodicity  and  coherence  individually, 
and  found  that  homomorphic  dereverberation  is  not  critically 
sensitive  in  either  case.   We  surmise  that  the  combination  of 
these  effects  would  not  be  catastrophic  to  performance.   Further 
investigation  is  merited  since  this  topic  has  an  important 
bearing  on  the  effectiveness  of  the  homomorphic  technique  in 
shallow  water  dereverberation  where  later  multiples  are 
significant. 

5.  Water  Travel  Time  Estimate 

In  practice,  passband  location  must  be  determined  by 
estimation  of  water  column  travel  time.   It  is  apparent  from 
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the  preceding  discussion  that  the  importance  of  accurate 
bottom  estimation  depends  on  the  seismic  environment.   In 
particular,  we  have  seen  that  MSR  and  bottom  characteristics  seen 
to  affect  cepstral  energy  distribution.   Proximity  of  multiples 
to  important  reflectors  also  dictates  resolution  constraints 
on  travel  time  estimation.   For  the  cases  analyzed  here  the 
required  travel  time  estimation  accuracy  would  be  20-40  msec 
since,  as  determined  in  the  foregoing  discussion,  a  stopband 
of  40-80  msec  is  generally  required  for  effective  dereverberation , 
Travel  time  estimation  error  of  more  than  half  the  stopband 
width  causes  the  major  multiple  contributions  to  be  excluded 
from  the  stopband.   This  degree  of  accuracy  can  easily  be 
attained  with  existing  bottom  tracking  algorithms. 

The  reflector/multiple  resolution  Implied  by  the  required 
stopband  widths  is  also  about  40-80  msec  for  the  data  tested. 
A  discussion  and  several  examples  of  reflector/multiple  resolu- 
tion are  included  in  the  following  section. 

6  Reflector/Multiple  Overlap 

The  following  six  figures  illustrate  homomorphic  dereverbera- 
tion of  signals  in  which  reflectors  and  multiples . are  closely 
situated. 

Figure  52a  shows  a  seismogram  with  reflectors  at  1.6,  2.02 
and  3.0  seconds.   The  first  and  second  multiples  directly  inter- 
fere with  reflectors  and  the  third  multiple  is  barely  visible 
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at  4.0  seconds.   The  result  of  applying  a  cepstral  stopband  at 
.96-1.02  seconds  is  seen  in  figure  52b.   Most  of  the  energy- 
near  two  seconds  has  been  removed  but  a  small  reflection  is 
still  visible.   The  signal  of  figure  52c  results  from  a  stop- 
band  of  .98-1.06  seconds  which  spans. the  time  domain  reflector 
location  as  well  as  the  multiple  location.   Residual  energy 
is  still  present  after  this  processing  which  indicates  that 
some  of  the  reflector  and  first  multiple  cepstrum  contributions 
are  distributed  beyond  the  immediate  vicinity  of  the  time 
domain  locations . 

The  seismogram  of  figure  53a  contains  reflectors  at  1.6, 
1.95  and  3.0  seconds.   Each  reflector  is  clearly  evident  after 
processing  (.98-1.06  second  stopband)  and  first  multiple  energy 
has  been  greatly  reduced  as  shown  in  figure  53b.   Some  of  the 
removed  energy  may  be  due  to  the  1.9  5  second  reflector;  however, 
the  first  multiple  in  this  example  is  very  coherent  and  its 
energy  is  more  likely  to  be  localized  near  1.0  second  in  the 
cepstrum.   Extension  of  the  stopband  closer  to  the  reflector 
location  (. 96-1. 06  second)  causes  complete  removal  of  the 
reflector  as  seen  in  figure  52c.   Visible  reduction  of  the 
second  multiple  at  3 . 0  seconds  and  an  internal  multiple  at  2 . 7 
seconds   is  apparent  in  figure  52b. 

The  seismogram  of  figure  54a  has  the  same  reflector  place- 
ment as  that  of  figure  53a  but  the  multiples  in  this  case  are 
smaller  and  less  coherent.   Application  of  a  cepstral  stopband 
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at  1.00-1.06  second  yields  Cfigure  54b)  visible  reduction 
of  all  three  multiples  and  well  resolved  reflections  at  the 
proper  locations . 

Figure  55a  illustrates  a  seismogram  with  considerable 
multiple  distortion  and  reflectors  at  1.5,  2.1  and  3.0  seconds. 
A  wide  cepstral  stopband,  .98-1.08  second,  has  a  very  minor 
effect  on  the  first  multiple  as  seen  in  figure  55b.   Extension 
of  the  stopband  to  1.12  seconds  causes  removal  of  virtually 
all  signal  energy  in  the  region.   The  relatively  wide  (.1  second) 
reflector/multiple  separation  cannot  be  effectively  exploited 
in  this  case  because  the  very  incoherent  multiple  in  a  low  MSR 
signal  requires  a  large  stopband  for  effective  removal  (see 
figure  4  8) . 

The  result  of  longpass  filtering  the  cepstrum  of  the  signal 
in  figure  55a  is  shown  in  figure  56.   All  cepstral  contributions 
prior  to  1.02  seconds  have  been  set  to  zero.   The  bottom  and 
later  reflectors  are  clearly  visible  but  the  1.5  second  reflector 
has  been  removed.   This  illustrates  one  drawback  to  longpass 
dereverberation  in  deeper  water.   This  effect  may  be  acceptable 
in  some  situations,  however.   For  instance,  good  resolution  of 
closely  spaced,  deep  reflectors  may  be  obtained  by  longpass 
filtering  whereas  the  earlier  reflectors  are  frequently  obvious 
before  processing. 

Figure  57  shows  longpass  results  for  a  signal  with  impulsive 
multiples  and  very  sharp  reflectors  at  2.2  and  3.0  seconds. 


Rl 


-123- 
R2 


Bl 


(a) 


B2,    R3 


B3 


WWyV~~ -wvw — /|/ia~^ 


(b) 


W^V 


(c) 


-vy\/w— 


1.0 


2.0 


3.0 


4.0 


5.0 
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Figure  56   (a)  Unprocessed  seismogram 

(b)  Result  of  longpass  filtering  cepstrum;  cutoff 
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The  stopband  limit  is  1.01  seconds  in  this  case.   Multiple 
removal  is  complete  and  both  reflectors  are  clearly  visible. 
Longpass  filtering  was  found  to  be  most  effective  for  signals 
of  this  type.   Very  sharp  reflectors  are  more  clearly  visible 
than  more  distorted  reflectors  of  comparable  energy,  especially 
in  noisy  signals.   Longpass  filtering  of  noisy  signals  is 
discussed  later  in  this  chapter. 

The  low  frequency  noise  present  in  the  reflector  region 
of  figure  57  is  typical  of  that  observed  in  several  longpass 
results.   The  cepstra  of  figures  49  and  50  contain  similar 
components.   No  attempt  was  made  in  this  analysis  to  remove 
this  type  of  noise.   The  reason  for  its  presence  has  not  been 
determined. 

7.   Additive  White  Noise  and  Filtering 

The  effects  of  additive  noise  are  not  addressed  in  the 
formulation  of  the  homomorphic  system  for  deconvolution  since 
there  are  no  apparent  cepstral  properties  of  noise  which  can  be 
exploited  for  its  removal.   Linear  filtering  is  a  more  suitable 
way  of  reducing  additive  noise  in  individual  signals.   This  can 
be  performed  prior  to,  or  in  conjunction  with,  homomorphic 
deconvolution.   The  effects  of  this  combined  processing  in  the 
special  case  of  Gaussian  white  noise  have  been  investigated  and 
are  discussed  here.   The  data  presented  represent  a  relatively 
small  number  of  experiments  performed  with  synthetic  data  and 
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noise  generated  as  previously  discussed.   These  results  cannot 
be  generalized  categorically  due  to  the  relatively  narrow 
scope  of  the  experiments  and  the  lack  of  precise  theoretical 
characterization  of  noise  properties  under  homomorphic  trans- 
formation.  The  results  do  indicate  performance  trends  and 
computational  effects  due  to  additive  noise  and  filtering. 

Since  noisy  signals  usually  undergo  linear  filtering  prior 
to  dereverberation,  we  begin  by  discussing  an  important 
computational  issue  which  arises  when  signals  are  bandpass 
filtered.   Such  filtering  introduces  spectral  regions  (stopbands) 
containing  little  or  no  spectral  energy.   Recall  that  the  homo- 
morphic transformation  involves  computing  the  logarithm  of  the 
Fourier  transform  of  the  signal.   Since  the  logarithm  of  zero 
is  not  defined,  this  computation  is  not  possible  in  frequency 
bands  where  the  Fourier  transform  is  zero.   In  the  case  of 
lowpass  filtered  signals  this  problem  can  frequently  be  overcome 
by  resampling  after  filtering.   Figures  58  a,  b  and  c  illustrate 
this  process.   If  the  sampling  frequency  is  decreased  to  the 
Nyquist  rate  implied  by  the  filter  cutoff  frequency,  the 
resulting  discrete  Fourier  transform  will-  include  only  the 
frequency  components  in  the  passband  region.   Figure  5  8d  shows 
a  distribution  of  spectral  energy  which  is  not  readily  amenable 
to  elimination  of  the  zero  region.   In  this  case  the  baseband 
region  is  zero  so  that  resampling  would  not  be  effective. 
Investigation  of  this  problem  is  currently  in  progress 
(Tribolet  [11] ) . 
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Figure  58   (a)  Unfiltered  spectrum  sampled  at  the  Nyquist  rate 

(b)  Spectrum  of  (a)  after  lowpass  filtering  at  u>c 
without  resampling. 

(c)  Filtered  spectrum  after  reampling. 

(d)  Spectrum  with  zeros  in  the  baseband  region. 
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It  was  found  in  this  analysis  that  computing  the  cepstrum 
of  an  over- sampled  signal  is  an  unstable  procedure  which,  at  best, 
requires  extensive  computation  time.   The  low  energy  regions 
of  the  spectrum  lead  to  spurious  phase  derivative  values. 
Hence,  the  phase  unwrapping  algorithm  proceeds  at  small  step 
sizes,  requiring   computation  of  many  intermediate  values  of  the 
discrete  Fourier  transform.   In  many  cases  the  computer  algorithm 
could  not  produce  an  unwrapped  phase  which  was  in  acceptable 
agreement  with  the  principal  value.   Cepstra  of  over-sampled 
signals  which  were  successfully  computed,  however,  generally 
yielded  dereverberation  results  comparable  to  those  of  properly 
sampled  signals  (see  table  2) . 

Addition  of  white  noise  was  found  to  have  a  computational 
effect  similar  to  that  described  above;  phase  unwrapping  time 
is  significantly  increased.   Heavy  weighting  (w=.98-.99) 
was  found  to  reduce  computation  time  considerably  for  noisy 
signals,  although  some  signals  require  small  phase  integration 
step  sizes  in  isolated  sections  even  when  substantial  weighting 
is  applied.   Recall  from  Chapter  II  that  there  is  no  assurance 
that  the  log-spectrum  is  adequately  sampled,  even  after  lowpass 
filtering  of  the  signal.   Hence,  the  integration  of  the  phase 
derivative  cannot  be  expected  to  proceed  quickly  in  all  cases. 

Smoothing  of  the  phase  derivative  was  attempted  to 
compensate  for  the  effects  of  noise.   A  three-point  moving 
average  was  applied  prior  to  integration.   The  resulting  inverse 
cepstra  bore  no  resemblance  to  the  original  seismograms. 
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FRACTION  OF  ENERGY  REMOVED 


SAMPLING 

INTERVAL   FILTER 
(MSEC)    CUTOFF (Hz)  1st  MULT   2nd  MULT   3rd  MULT 


4.88 

50 

.91 

.064 

.25 

9.77 

50 

-.02 

.12 

.24 

19.53 

50 

.94       .044 

.25 

4.88 

20 

.94 

.083 

.15 

9.77 

20 

.93       .48 

.13 

19.53 

20 

.98       .12 

.55 

.  .. 

■ 

TABLE  2 


EFFECT  OF  RESAMPLING  ON  MULTIPLE  REMOVAL  FOR  A 
NOISELESS  SEISMOGRAM  LOWPASS  FILTERED  AT  5  0  and 
20  Hz.   DE REVERBERATION  WAS  ACCOMPLISHED  BY  APPLY- 
ING AN  80  MSEC  CEPSTRAL  STOPBAND  AT  THE  FIRST 
MULTIPLE  LOCATION. 
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Figure  59  summarizes  the  observed  effects  of  noise  on 
first  multiple  removal  by  homomorphic  processing.   Each 
seismogram  was  lowpass  filtered  at  50  Hz  and  resampled  by  a 
factor  of  four  prior  to  multiple  removal.   Percentage  of 
multiple  energy  removed  shows  a  consistent  decrease  with  increas- 
ing noise  level.   The  rates  of  decrease  and  amounts  of  energy 
removed  are  seen  to  vary  widely  from  one  signal  to  another. 
Examples  of  noisy  signals  before  and  after  processing  are  shown 
in  figures  60  and  61.   Significant  reduction  of  the  first  multiple 
is  evident  in  both  examples.   In  the  first  case,  the  noise  level 
is  moderate  and  multiple  reduction  has  a  marked  effect  on  visual 
quality  of  the  signal.   The  noise  level  in  figure  61  is  consider- 
ably higher,  resulting  in  marginal  improvement  due  to  dereverbera- 
tion. 

Second  and  third  multiple  energy  removed  was  found  to 
decrease  generally  with  decreasing  noise  also,  and  in  some 
high  noise  cases  the  later  multiples  were  actually  enhanced. 
Data  for  a  typical  signal  are  tabulated  below. 


ENERGY  REMOVED 


STANDARD 

DEVIATION 

OF  NOISE    2nd  MULT.  3rd  MULT 


0 

.29                    .32 

25 

.19 

.15 

50 

.04 

.02 

100 

-.21 

-.17 

TABLE  3 
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Figure  59   Effect  of  added  noise  on  homomorphic  dereverberation 
All  signals  are  lowpass  filtered  at  50  Hz  and 
resampled. 
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Figure  60   (a)  Seismogram  with  moderately  high  noise,  F  =50  Hz. 
(b)  Result  of  cepstral  notch  filtering  (.8-154  sec). 
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Figure  61   (a)  Sei sinogram  with  very  high  noise,  F  =50  Hz. 

(b)  Result  of  cepstral  notch  filtering0 (. 80-1 . 04  sec) 
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Reference  to  figure  61  reveals  that  for  the  higher  noise 
amplitudes  the  later  multiples  are  buried  so  that  their  removal 
is  not  important. 

Reflector  distortion,  summarized  for  a  typical  case  in 
table  4,  was  found  to  be  generally  more  severe  than  in  noise- 
less signals  but  rarely  more  than  25%. 


STANDARD 
DEVIATION 
OF  NOISE 


.REFLECTOR  ENERGY  REMOVED 
12  3 


0 

-.18 

.02 

.06 

25 

-.22 

•     -.25 

-.22 

50 

-.07 

-.14 

-.05 

100 

.14 

-.12                -.09 

• 

TABLE  4 

EFFECT  OF  NOISE  ON  REFLECTOR  DISTORTION 
FOR  A  TYPICAL  SEISMOGRAM  WITH  3  REFLECTORS, 
FILTERED  AT  50  Hz. 


The  above  results  indicate  that  the  addition  of  noise 
leads  to  decreasing  performance  with  respect  to  all  three 
criteria.   The  behavior  is  somewhat  spurious  and  frequently 
exhibits  wide  variation  from  signal  to  signal. 
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Several  lowpass  filter  bandwidths  (20,  30,  50,  70,  and 
90  Hz)  were  applied  to  the  signals  evaluated  in  figure  59. 
The  results  are  shown  in  figure  62.   In  each  case,  first  multiple 
removal  is  least  effective  when  the  signals  are  prefiltered  at 
50  Hz.   Figures  63  and  64  illustrate  this  behavior.   The  same 
seismogram  is  shown  in  figures  6  3a  and  64a,  lowpass  filtered 
at  50  and  20  Hz  respectively.   The  difference  in  appearance  is 
dramatic.   Resolution  of  the  reflectors  on  either  side  of  the 
first  multiple  (at  1.55  and  2.4  sec)  is  greatly  improved  by 
dereverberation  in  the  latter  case  while  figure  6  3  shows  very 
little  visual  improvement. 

The  behavior  illustrated  in  these  figures  cannot  be  fully 
explained  on  the  basis  of  the  data  available.   The  signals 
tested  have  dissimilar  reflector  locations  and  varying  amounts 
of  multiple  distortion  which  implies  that  the  similar  performance 
dips  are  not  due  to  similarities  in  signal  configuration.   The 
observed  behavior  may  be  due  to  the  properties  of  the  source 
signature  (which  is  common  to  all  three  seismograms)  or  the 
characteristics  of  the  recursive  third  order  Butterworth 
filter  employed.   More  data  using  different  filter  routines 
and  a  wide  range  of  signal  characteristics  is  needed  to 
completely  characterize  this  behavior.   The  results  obtained 
here  suggest  that  a  filter  bandwidth  very  close  to  the  bandwidth 
of  the  noiseless  signal  leads  to  the  best  dereverberation 
performance.   The  same  result  was  obtained  for  the  TDL  algorithm. 
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Figure  62   Effects  of  lov/pass  filtering  signals  with  added 
white  noise  prior  to  homomorphic  dereverberation 
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Figure  63   (a)  Noisy  seismogram  lowpass  filtered  at  50  Hz.  and 

resampled  at  51  Kz. 
Cb)  Result  of  cepstral  notch  filtering;  stopband 
.8-1.04  sec. 
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Figure  64   Noisy  seism.ogram  lowpass  filtered  at  20  Hz  and 

resampled  at  51  Hz.   (b)  Result  of  cepstral  filtering; 
stopband  .8-1.04  sec. 
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Additive  noise  was  found  to  have  a  particularly  adverse 
effect  on  longpass  filtered  signals.   Noise,  which  may  be 
relatively  low  in  the  received  signal,  tends  to  be  amplified 
with  respect  to  reflectors  in  the  process  of  longpass  filtering, 
especially  in  the  later  portion  of  the  signal.   Although  long- 
pass  filtering  removes  a  large  part  of  the  noise  energy  with 
the  source,  the  overall  effect  is  usually  a  decrease  in  SNR. 
The  advantage  of  longpass  filtering,  which  includes  source 
deconvolution,  is  that  greater  resolution  of  close  reflectors 
can  be  achieved.   It  was  found  that  this  approach  is  worthwhile 
in  low  noise  signals  but  not  effective  when  the  white  noise 
level  is  significant  with  respect  to  reflector  amplitudes. 

Proper  resampling  after  prefiltering  was  found  to  be  very 
important  for  successful  longpass  filtering.   Figure  65a 
illustrates  the  effect  of  longpass  filtering  the  cepstrum  of  a 
noisy  signal  which  was  first  lowpass  filtered  at  20  Hz  but 
not  resampled.   The  filtered  signal  contains  relatively  low 
noise  and  the  sampling  frequency  is  205  Hz.   The  processed 
result  of  figure  65b  is  useless  due  to  the  high  sampling  rate. 
Reduction  of  the  sampling  frequency  to  102.5  Hz  leads  to  the 
processed  signal  of  figure  66a.   The  3.5  second  reflector  is 
clear  but  the  high  background  noise  almost  obscures  the  2.7 
second  reflector.   Multiple  removal  is  complete.   Resampling 
to  51  Hz,  which  is  approximately  the  Nyquist  rate  in  this  case, 
leads  to  some  improvement  (figure  66b)  but  the  SNR  is  much  lower 
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Figure  65   (a)  Noisy  seismogram  lowpass  filtered  at  20  Hz  but 

not  resampled. 
(b)  Result  of  longpass  filtering  the  cepstrura  of  (a) 


-142- 


Rl 


R2 


81 


B2 


93 


■vs**^ — Vir^i^^^ 


(a) 


^^^Vv\A^^ 


(b) 


1.0 


2.0 


3.0 


4.0 


5.0 


Figure  66   Longpass  processing  results  (a)  for  the  signal  of 

figure  65a,  resampled  at  102.5  Hz  before  processing. 
(b)  for  the  signal  of  figure  65a,  re sampled  at  51  Hz 
before  processing  (c)  for  a  signal  identical  to 
figure  65a  with  one  half  the  white  noise  level, 
resampled  at  51  Hz. 
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than  before  dereverberation  (figure  65a) ,   Figure  66c  results 
from  processing  identical  to  that  of  figure  6  6b,  on  the  same 
signal,  with  the  noise  amplitude  halved.   Both  reflectors  are 
clear  and  the  multiples  have  been  removed  but  the  background 
noise  is  much  higher  even  than  in  the  signal  of  figure  65a. 
Thus,  longpass  filtering  of  noisy  signals  is  seen  to  involve 
a  trade-off  between  effective  dereverberation  and  decrease  in 
SNR.   For  signals  with  moderate  to  heavy  noise  the  reduction  in 
SNR  was  found  to  be  unacceptable  in  the  cases  tested. 

Although  this  subject  was  not  extensively  investigated 
there  is  some  indication  that  lowpass  filtering  can  be  employed 
to  improve  the  results  of  subsequent  longpass  processing. 
Figures  67  and  68  illustrate  the  longpass  processing  of  a 
noisy  signal  after  lowpass  prefiltering  at  (a)  70  Hz,  (b)  50  Hz 
and  (c)  30  Hz.   In  each  of  the  signals  in  figure  68  the  multiples 
have  been  removed  but  the  reflector  resolution  improves 
considerably  from  (a)  to  (c) .   These  results  are  reasonable 
in  that  improvement  of  performance  coincides  with  increasing 
rejection  of  out-of-band  noise;  however,  the  filter  and  signal 
characteristics  must  be  studied  more  closely  to  explain  this 
behavior  accurately. 

D.   Comparative  Examples  of  Processing  Results 

The  foregoing  results  and  discussion  illustrate  the 
performance  of  the  TDL  and  homomorphic  dereverberation  algorithms 
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Figure  67   Noisy  seismograrn  lowpass  filtered  at  three  different 
frequencies  (a)  70  Hz  (b)  50  Hz  (c)  30  Hz.   Each  has 
been  re samp led  at  51  Hz  which  is  the  approximate 
Nyquist  rate  for  the  signal  without  noise. 
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Figure  68   Results  of  longpass  filtering  the  seismograras  of 
figure  65  a,  b  and  c,  respectively. 
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for  a  variety  of  processing  conditions.   Several  trends  in 
performance  are  apparent.   Before  proceeding  to  a  summary 
and  discussion  of  the  relative  strengths  and  weaknesses,  we 
present  several  examples  which  allow  direct  visual  comparison 
of  the  techniques.   Each  of  the  following  figures  contains 
(a)  an  unprocessed  seismogram  and  the  two  processed  results 
obtained  by  (b)  homomorphic  and  (c)  TDL  filtering. 

Figure  69  illustrates  a  noiseless  seismogram  with  multiple/ 
reflector  overlap  at  both  2  and  3  seconds.   Both  algorithms 
leave  a  small  amount  of  energy  at  the  first  location  and  effect 
only  slight  reduction  at  the  second.   In  general,  both  methods 
were  found  to  eradicate  reflectors  which  are  extremely  close 
to  the  first  multiple  and  retain  signal  components  which  closely 
coincide  with  later  multiples. 

In  figure  70  the  multiple  onset  occurs  .1  second  before 
the  reflector  and  the  MSR  is  considerably  lower  than  in  the 
previous  figure.   In  this  case  the  separation  is  great  enough 
that  both  methods  retain  a  significant  amount  of  signal  energy 
near  2.1  seconds.   The  homomorphic  result  is  considerably 
sharper  although  very  little  energy  has  been  removed.   Multiple/ 
reflector  separation  of  .1  second  was  found  to  be  the  approxi- 
mate resolution  limit  of  both  techniques  when  reflector  onset 
is  later  than  multiple  onset  and  travel  time  is  estimated 
accurately. 

Figure  71  shows  a  reflector  at  1.95  seconds,  slightly 
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Figure  69   (a)  Unprocessed  seismogram  with  reflectors  at  1.5,  2.02 

and  3.0  seconds. 

(b)  Result  of  homomorphic  processing. 

(c)  Result  of  TDL  processing. 
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Figure  70   (a)  Unprocessed  seismogram  with  reflectors  at  1.5,  2.1 

and  3.0  sec. 

(b)  Result  of  homomorphic  processing 

(c)  Result  of  TDL  processing. 
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Figure  71   (a)  Unprocessed  seismogram  with  reflectors  at  1.6,  1.95 

and  3.0  sec. 

(b)  Result  of  homomorphic  processing. 

(c)  Result  of  TDL  processing. 
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before  the  multiple.   Effective  dereverberation  is  accomplished 
in  both  b  and  c.   As  in  the  previous  figure  the  reflector  is 
better  resolved  by  homomorphic  processing.   Resolution  of  both 
techniques  was  generally  observed  to  be  slightly  better  when 
reflector  onset  is  earlier  than  multiple  onset,  provided  travel 
time  is  estimated  accurately.   In  such  cases  of  severe  inter- 
ference the  homomorphic  filtering  generally  yields  more  distinct 
reflections.   A  further  example  of  this  behavior  is  shown  in 
figure  72. 

The  following  three  figures  illustrate  dereverberation  of 
noisy  signals.   Each  seismogram  has  been  lowpass  filtered  at 
30  Hz  and  the  homomorphic  outputs  as  shown  have  been  resampled 
at  51  Hz.   In  this  first  case  (figure  73)  the  performance  of 
both  methods  is  comparable.   The  first  multiple  is  almost 
completely  removed  while  other  regions  of  the  signal  are  not 
visibly  affected.   Figure  74  also  shows  comparable  multiple 
removal,  however,  the  resolution  of  the  reflectors  at  2.6  and 
3.4  seconds  is  somewhat  better  after  TDL  filtering.   The 
homomorphic  algorithm  was  found  to  produce  higher  random  noise 
spikes  than  the  TDL  filter.   This  effect  is  present  in  figure 
74  and  again  in  figure  75.   In  both  examples  the  homomorphic 
method  achieves  slightly  better  multiple  removal  but  the 
overall  noise  level  in  the  result  appears  higher. 
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Figure  72   (a)  Unprocessed  seismogram  with  reflectors  at  1.6, 

1.95  and  3.0  sec. 
Cb)  Result  of  homomorphic  processing, 
(c)  Result  of  TDL  processing. 
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Figure  73   (a)  Noisy  seismogram  with  reflectors  at  2.7  and 

3.5  sec,  and  lowpass  filtered  at  30  Hz. 

(b)  Result  of  homomorphic  processing. 

(c)  Result  of  TDL  processing. 
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Figure  74   Ca)  Noisy  seismogram  v;ith  reflectors  at  1.55,  2.4  and 

3.25  sec,  lowpass  filtered  at  30  Hz. 

(b)  Result  of  homomorphic  processing. 

(c)  Result  of  TDL  processing. 
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Figure  75   Ca)  Noisy  seismogran  with  reflectors  at  2.6  and  3.4 

sec,  lovTpass  filtered  at  3  0  Hz. 

(b)  Result  of  homomorphic  processing. 

(c)  Result  of  TDL  processing. 
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CHAPTER   V 
DISCUSSION    AMD    SUMMARY 

In  the  preceding  chapters  we  have  (1)  examined  the 
theoretical  structures  of  the  TDL  and  homomorphic  dereverbera- 
tion  techniques ,  (2)  established  a  comprehensive  set  of  perfor- 
mance criteria,  and  (3)  presented  the  results  of  applying  both 
methods  to  synthetic  data.   Our  approach  has  been  essentially 
that  of  perturbation  analysis.   Through  variation  of  environ- 
mental and  signal  processing  parameters  we  have  observed 
performance  trends  due  to  deviations  from  the  ideal  theoretical 
models  upon  which  the  methods  are  based.   In  a  more  practical 
sense  the  parameter  variations  simulate  a  range  of  seismic 
processing  conditions.   Since  there  are  many  different  environ- 
ments in  which  these  algorithms  may  be  applied,  we  have  not 
emphasized  the  absolute  performance  figures  obtained  from  the 
simple,  synthetic  data  utilized  here.   Rather,  we  have  tried 
to  present  a  behavior  profile  of  both  algorithms  which  indicates 
the  basic  trends  and  sensitivities  with  respect  to  a  number 
of  parameters,  interpreted  in  terms  of  their  theoretical 
structures  and  assumptions.   This  approach  is  intended  to  give 
a  more  general  indication  of  the  dereverberation  performance 
to  be  expected  in  different  situations.   We  conclude  with  a 
summary  and  discussion  of  the  comparative  results  presented  in 
Chapter  IV. 
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In  low  noise,  deep  water  seismograms  the  methods  have 
been  found  to  be   comparable  for  reducing  dominant  first 
multiples  (by  75-95%  in  most  cases) .   The  TDL  filter  requires 
accurate  signal  statistics  including  an  estimate  of  the  multiple- 
reflector  crosscorrelation  function,  which  must  be  approximated. 
This  leads  to  degraded  performance  in  shallow  water  situations 
where  significant  reflector  energy  is  within  the  crosscorrelation 
window.   The  homomorphic  method  requires  no  statistical  character- 
ization of  the  signal  and  thus  has  no  similar  performance 
degradation  in  shallow  water;  however,  the  three-stage  cepstral 
transformation  requires  extensive  computation  which  may  be  an 
important  limitation  for  at-sea  processing  systems.   (This 
issue  will  be  discussed  in  more  'detail  later.) 

Although  the  effects  of  aperiodicity  could  not  be  thoroughly 
evaluated  experimentally  the  derived  result  expressed  in 
equation  (11)  suggests  that  the  homomorphic  algorithm  has  the 
potential  to  reduce  later,  aperiodic  multiples.   The  combina- 
tion of  such  processing  with  source  deconvolution  by  longpass 
filtering  the  cepstra  of  shallow  water  seismograms  appears  to 
be  the  most  promising  application  of  homomorphic  dereverbera- 
tion.   More  extensive  practical  evaluation  is  needed  in  this 
area. 

Closely  spaced,  aperiodic  multiples  destroy  the  coherence 
of  the  approximated  crosscorrelation  function  at  shifts  near 
the  two-way  travel  time  which,  again,  limits  the  effectiveness 
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of  the  TDL  filter  in  shallow  water  signals. 

Increasing  white  noise  level  causes  monotonically  decreas- 
ing performance  in  both  methods  as  illsustrated  in  figures  38 
and  57.   The  TDL  performance  fell  off  more  slowly  with  noise 
when  both  techniques  were  applied  to. similar  signals.   Bandpass 
filtering  leads  to  a  consistently  higher  percentage  of  multiple 
reduction  by  TDL  processing  of  noisy  signals.   Filter  effects 
on  homomorphic  processing  are  more  complex.   The  cases  evaluated 
indicate  that  multiple  energy  removed  is  not  a  monotonic  func- 
tion of  filter  cutoff  frequency.   Considerably  more  data  are 
required  to  determine  the  precise  effects  of  filter  bandwidth 
and  phase  characteristics.   The  resampling  which  was  found  to 
be  helpful  after  bandpass  filtering  may  have  a  .detrimental 
effect  on  visual  record  quality,  so  that  interpolation  may  be 
desirable  in  some  cases. 

Reflector  distortion  does  not  appear  to  be  a  problem  for 
either  technique  except  in  cases  where  overlap  is  severe.   In 
most  cases  tested  less  than  10%  of  the  reflector  energy  was 
removed.   Close  proximity  of  reflectors  to  the  first  multiple 
frequently  leads  to  severe  distortion  by. both  methods  due  to 
the  resulting  bias  of  the  crosscorrelation  function  and  the 
lack  of  sufficient  cepstral  separation.   Reflectors  close  to 
later  multiples  are  usually  well  preserved.   This  behavior 
suggests  two  applications  of  the  homomorphic  algorithm.   First, 
when  regions  of  geological  interest  occur  after  the  onset  of 
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the  first  multiple,  a  wide  cepstral  stopband  can  be  employed  at 
the  first  multiple  location  to  reduce  later  multiples  without 
distorting  reflectors.   In  very  shallow  water  the  stopband  may 
be  extended  to  include  more  than  one  multiple.   A  second 
possibility  for  avoiding  reflector  distortion  is  the  use  of 
weighting  coefficients  greater  than  1.0  to  exploit  the  properties 
of  mixed  phase  sequences.   The  object  of  this  weighting  is  to 
make  the  reflector  train  mixed  phase  while  keeping  the  multiple 
sequence  minimum  phase.   This  appears  to  be  feasible  in  many 
situations  since  the  z-plane  zero  of  the  multiple  sequence  is 
usually  well  inside  the  unit  circle.   Moving  some  of  the  reflector 
train  zeros  outside  of  the  unit  circle  (i.e.,  making  it  mixed 
phase)  will,  in  general,  cause  some  of  the  cepstral  energy  due 
to  the  reflectors  to  occupy  the  negative  quefrency  region. 
Even  if  the  reflector  train  has  a  maximum  phase  component 
before  weighting  the  same  effect  can  be  expected.   Thus,  the 
amount  of  reflector  energy  near  the  first  multiple  location 
may  be  reduced.   Although  there  is  no  guarantee  that  the 
resulting  notch  filtered  cepstrum  will  transform  to  a  seismogram 
with  less  distortion,  this  technique  appears  to  be  worthy  of 
investigation . 

We  recall  one  other  reflector  distortion  effect  which  was 
seen  in  Chapter  IV.   We  saw  in  figure  55  that  dereverberation 
by  longpass  filtering  completely  removes  reflectors  occurring 
prior  to  multiple  onset.   It  was  noted  that  this  effect  may  be 
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acceptable  if  it  leads  to  better  resolution  of  smaller,  deep 
reflectors . 

In  terms  of  our  third  criterion,  visual  improvement  of  the 
signal,  both  techniques  were  seen  to  have  advantages  and 
disadvantages.   Homomorphic  processing  usually  resulted  in 
better  resolution  of  interfering   signal  components  in  the  low 
noise  signals  processed.   Increasing  noise,  however,  was  seen 
to  cause  homomorphic  results   more  prone  to  random  noise  spikes 
which  degrade  the  interpretability  of  the  record.   The  results 
of  TDL  filtering  had  somewhat  better  visual  resolution  in  the 
noisier  signals  processed.   As  noted  previously  the  visual 
advantage  of  the  TDL  in  this  case  is  partially  due  to  the 
noisier  appearance  of  the  resampled  signals  produced  by 
homomorphic  processing. 

Longpass  filtering  was  seen  to  provide  the  most  effective 
dereverberation,  the  best  reflector  resolution  and  the  best 
overall  visual  quality  in  ideal  cases.   Unfortunately  it 
degenerates  quickly  with  noise  and  could  not  be  successfully 
applied  to  very  noisy  signals  or  signals  with  important  geological 
regions  preceding  the  multiple  onset.   Further  research  and 
experience  may  well  lead  to  more  extensive  applicability  of 
this  technique. 

The  relative  simplicity  of  the  TDL  algorithm  makes  it 
much  more  desirable  from  a  computational  standpoint.   TDL 
dereverberation  of  a  1024  point  signal  can  be  accomplished  in 
3  seconds  or  less  for  the  operator  lengths  typically  required. 
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The  three  major  computational  steps  are  the  correlation 
operations,  solution  of  the  estimator  equations  and  convolu- 
tion of  the  operator  with  the  signal. 

The  homomorphic  computations  include  weighting,  four  FFT's 
computation  of  the  complex  logarithm,  phase  unwrapping,  linear 
filtering,  complex  exponentiation  and  unweighting.   This 
algorithm  can  be  expected  to  take  2  0  seconds  or  more  for  a 
1024  point  sequence  on  a  small  processing  computer.   In  this 
analysis  the  phase  unwrapping  computation  took  over  one  minute 
for  some  noisy  signals.   These  figures  are  highly  dependent 
on  hardware  available  and  programming  efficiency  but,  in 
general,  homomorphic  dereverberation  is  several  times  slower 
than  the  TDL  algorithm.   Special  purpose  hardware  could  be 
used  to  reduce  homomorphic  computation  tine  significantly, 
but  the  method  has  not  been  implemented  for  processing  on  a 
large  scale  thusfar. 

Storage  requirements  for  the  homomorphic  algorithm  vary 
with  the  FFT  routine  used,  method  of  cepstrum  computation  and 
cepstrum  length.   The  program  used  for  this  analysis  requires 
about  12  *  N  bytes  of  core  and  4  *  N  bytes  of  disc  storage, 
where  N  is  the  cepstrum  length.   N  was  twice  the  signal  length 
in  the  cepstra  computed.   Shorter  cepstrum  lengths,  as 
determined  by  trial-and-error ,  may  produce  results  with 
acceptably  low  aliasing  in  many  cases.   The  TDL  dereverberation 
program  used  requires  about  6  *  L  bytes  of  core,  where  L  is 
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the  data  sequence  length.   No  disc  storage  is  required  in 
this  computation.   These  stbrage  requirements  apply  to  a  float- 
ing point  processing  scheme  on  a  machine  (HP-2100)  which  uses 
four  byte  floating  point  words. 

In  conclusion,  we  make  some  general  observations  concern- 
ing the  results  of  this  analysis. 

The  TDL  dereverberation  scheme  is  a  simple  and  efficient 
technique  which  has  demonstrated  effectiveness  in  removing 
deep  water  multiples.   The  analytical  structure  is  well  under- 
stood and  its  performance  characteristics  have  been  explained 
here  in  terms  of  that  structure.   Further  refinements  in  its 
implementation  may  be  possible  but  its  potential  is  essentially 
clear  at  this  point. 

Homomorphic  dereverberation  is  complex,  relatively  untested 
and  requires  extensive  computation.   It  has  been  shown  here 
to  be  effective  on  synthetic  data.   It  appears  to  be  particularly 
promising  for  shallow  water  dereverberation.   The  complexity 
of  the  method  leads  to  a  number  of  possible  analytic  and 
computational  techniques  which  can  be  utilized  in  its  applica- 
tion.  Certainly,  its  full  potential  has  not  yet  been  determined 
and  further  investigation  is  warranted. 
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APPENDIX   A 
COMPUTATION    OF    THE    PHASE    DERIVATIVE    OF    THE    Z-TRANSFORM 

When  computing  the  complex  cepstrum  of  a  sequence,  x(n), 
it  is  necessary  to  determine  the  unique,  continuous  phase  of 
X(z).   One  way  of  obtaining  the  continuous  phase  is  to  first 
compute  its  derivative  and  then  integrate  numerically.   The 
computation  of  the  phase  derivative  from  x(n)  is  discussed 
in  detail  here. 

We  begin  with  the  z-transform  of  x(n), 

00 

X(z)  =  I      x(n)  z"n  =  X  (z)  +  jXjCz). 

n=-oo 

Taking  the  complex  natural  logarithm 

log  X(z)  =  X(z)  =  log  |x(z)|  +  j  arg  X(z) 

We  see  that  the  phase  of  X(z)  is  equal  to  the  imaginary 
part  of  its  natural  logarithm.   The  derivative  of  X-^z)  can  be 
expressed  in  terms  of  easily  computable  quantities. 

f*xU)  =  _d  log[x(z)]  =  2LLLZ)  .  (1) 

dz      dz  X(z) 

where  the  prime  indicates  differentiation  with  respect  to  z. 
Expanding  CD  and  solving  for  X-^z), 
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X'(z)  +  jX  (z) 

XR(z)  +  jX  (z)  =   _* 1 

XR(z)  +  jXjCz) 


X'(z)  -  jx'(z) 

XT(z)  =   -I ^—   +  j  Xp(z) 

X  (z)  +  jX  (z) 


Separating  the  RHS  into  real  and  imaginary  parts, 


X  (z)  X  (z)  -  X_(z)  Xp(z) 
XI(z)  =  -5 I 1 B_ 

XR2(z)  +  X*(z) 


+  j 


(XR(z)  X'(z)  +  XT(z)  X*  (z)) 
X_(z) ^ 5 1 1 

2        2 
XR(z)  +  Xj(z) 


The  real  part  yields  an  expression  for  the  phase  derivative, 


X_(z)  Xl(z)  -  XT(z)  X_(z) 
Xx(z)    =   -5 1 1 *_ 

XR(z)  +  X*(z) 


(2) 


Since  the  z-transforra  is  actually  evaluated  on  the  unit  circle 
using  the  discrete  Fourier  transform  (DFT)  we  set  z  =  e-^. 


X'x{e^)    =  -* 


93»)  X;(e^)  -  Xl(e^)  XR(e^) 


£(ejw)  +  X?(e^) 


(3) 


R 


Derivatives  with  respect  to  e^   may  be  replaced  by  -=—   since 
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dX(e3b)) 


dco 


=  ]  e 


30J 


dx(e^) 


de 


ju 


and  we  thus  have  a  common  factor  of  j  e^w  in  each  term  of  (3) . 

Hence,  we  need  only  compute  the  real  and  imaginary  parts 
of  X(e-^  )  and   -3—  (xCe-1^))  and  combine  them  as  indicated  in  (3) 
The  derivative  of  Xte-^)  is  easily  obtained  from  the  sequence 
n  x(n)  as  f ollows : 


v  ,  iw,    r      /x   -icon    .  ,  X(eJ  ) 
X(eJ  )  =  }   n  x(n)  e  J    =    jd 


n 


n= 


JU), 

doo 


Re[Xn(e:|0))]  = 


-m*3a) 


Im[Xn(e:|w)]  = 


X^(e^) 


The  required  computation  in  terms  of  the  DFT  is  then 


Xj(k) 


-(xR(k)  xnR(k)  +  x1(k)   xnI(k)) 


XR(k)  +  Xj(k) 
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